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We present a complete numerical study of cosmological models with a time dependent 
coupling between the dark energy component driving the present accelerated expansion of the 
Universe and the Cold Dark Matter (CDM) fluid. Depending on the functional form of the 
coupling strength, these models show a range of possible intermediate behaviors between the 
standard ACDM background evolution and the widely studied case of interacting dark energy 
models with a constant coupling. These different background evolutions play a crucial role in 
the growth of cosmic structures, and determine strikingly different effects of the coupling on 
the internal dynamics of nonlinear objects. By means of a suitable modification of the cosmo- 
logical N-body code GADGET-2 we have performed a series of high-resolution N-body simu- 
lations of structure formation in the context of interacting dark energy models with variable 
couplings. Depending on the type of background evolution, the halo density profiles are found 
to be either less or more concentrated with respect to ACDM, contrarily to what happens for 
constant coupling models where concentrations can only decrease. However, for some spe- 
cific choice of the interaction function the reduction of halo concentrations can be larger than 
in constant coupling scenarios. We also find that different types of coupling evolution deter- 
mine specific features in the growth of large scale structures, like peculiar distortions of the 
matter power spectrum shape or different time evolutions of the halo mass function. Further- 
more, also for time dependent couplings baryons and CDM develop a bias already on large 
scales, which is progressively enhanced for smaller and smaller scales, and the effect can be 
significantly larger compared to constant coupling scenarios. The same happens to the baryon 
fraction of halos, which can be more significantly reduced below its universal value in variable 
coupling models with respect to constant coupling cosmologies. In general, we find that time 
dependent interactions between dark energy and CDM can in some cases determine stronger 
effects on structure formation as compared to the constant coupling case, with a significantly 
weaker impact on the background evolution of the Universe, and might therefore provide a 
more viable possibility to alleviate the tensions between observations and the ACDM model 
on small scales than the constant coupling scenario. 
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1 INTRODUCTION 

According to the present interpretation of the vast amount of cos- 
mological data that have been collected over the last years, the Uni- 
verse in which we live has a nearly flat spatial geometry, it ex- 
pands with a rate of ~ 70 km s~^ Mpc~^, and the total amount 
of matter it contains accounts only for ~ 27% of the critical en- 
ergy density that is required to justify its spatial flatness. Further- 
more, the expansion seems to have entered an accelerated phase 
since about 6 billion years, and the missing ~ 73% of the critical 
energy density is assumed to be in the form of some dark energy 
(DE) component able to drive such accelerated expansion. This de- 



tailed picture can be obtained by combining a wide variety of differ- 
ent and complementary cosmological datasets, ranging fro m Cos - 
mic Microwave Backround (CM B) jKoinat su et al. 2009 , 1201 Ol), 
to La r ge Scale S tructure surveys jPercival et al.ll2001k ICole et al.l 
120051 ; iReid et al.i 2010). to observations of Type la Supemovae 



(Snia) JRiess et al.lll998l ; IPerlmutter et"ai]| 19991; lAstier etal .1120061; 
iKowals ki et al.' '2008*) and Baryon Acoustic Oscillations (BAO) 
(e.g. [Percival et al. 2010). While the matter fraction of the Universe 
is known to be composed only for ~ 15% by baryonic particles, 
with the remaining mass in the form of some collisionless, weakly 
interacting Cold Dark Matter (CDM) component, the nature of the 
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DE fraction is yet completely unknown, and its understanding con- 
stitutes one of the major challenges in modern cosmology. 

The simplest possibility of a cosmological constant A, whose 
energy density remains constant throughout the whole expansion 
history of the Universe, is in good agreement with a very large 
amount of cosmological and astrophysical data, and this is the rea- 
son for the establishment of the ACDM scenario as the present 
standard cosmological model. Nevertheless, the nature of the cos- 
mological constant raises two fundamental questions concerning 
the very finely tuned value of its energy density (the "fine tuning 
problem") and the beginning of its domination over CDM only at 
relatively recent cosmological epochs (the "coincidence problem"). 
For this reason, alternative possibilities of dynamically evolving 
DE components have been proposed, in particular models where 
the DE is ide ntified with a classical scalar field as fo r the case of 
quintessence ([W etterich 1988; Ratra & Peebleslll988h or k-essence 
( lArmendariz?icon et al. 2000. 200 10 . 

As a further extension of these dynamical DE scenarios, dif- 
ferent possible forms of interaction between the DE component 
and the matter sector of the Universe have been suggested and in- 
vestigated injhe_jiteraturejas_e^the_gene^^ 



abundance of satellites in CDM halos JNavarro et al.l Il996p. to 
the observed low baryon fraction in large galaxy clusters i Ettoril 



(see e .g. iKamenshch ik et al.ll 200ll; Isilic et al.| l2002l: IS ento et al.l 



I2OO2I; ICarturan & F inelli 20031: lAmendola et al. 2003a), unified 



dark matter models aMainini & Bonomettd 12004: Bertacca et al 
2007, 2010'), extend ed q uintessence models I Perrotta et al. 2000|; 
Baccigalupi et al. 2000; Pettorino e t al.l uOOSh, and coupled DE 
(,Wetterichil995. : ,Amendolai200a ,2004 IPettorino & Baccigalupi 
uOOSh . Any form of interaction between the DE sector and 
other matter species, like CDM or massive neutrinos (as in e.g. 
lAmendola et alj2008l) would leave distinctive features in the back- 
ground expansion history o f the Universe an d in the growth of 
cosmic structures (see e.g. iBrax et al.l uOlOl) which could pro- 
vide new ways to tackle the problem of the nature of DE. It 
is therefore essential to understand the impact that the interac- 
tion woul d have on observable quantities as e.g. the properties 
of CMB JAmendola erai]|2003bl : Umendola & Ouercellinill2003l : 



1 lo 



iBean et al. 2008; La Vacca et al. 2009^, of larg e scale structure for- 
matio n (Koivisto 2005; Bean et al. 2008; Bal di & Pettorincll2010l ; 
iBaldi & Viell 1201 Ol). and of the nonlinear newtonian dynamics 
at small scales jPerrottaet alJl2()04l;lMainini & Bonomettoll2006l ; 



iMaccio et aLll2004ISaracco et alJl2010l ; lBaldi et alj2010l) . 

In the present work, in particular, we consider the case of a 
cosmological scenario where the DE scalar field interacts with the 
CDM fluid with a coupling strength that evolves in time, therefore 
generalizing the very widely studied case of constant couplings for 
the DE interaction. Other forms of effectively time dependent cou- 
plings, where the interaction depends on linear or nonlinear com- 
binations of the energy densities of the interac ting fluids and of 
their t i me derivatives have b e en stu d eid in e.g.lBarrow & Cliftonl 
( I2OO6I) ; ICaldera-Cabral et alj ( |2009|) ; IChimentd 1 I2OIOI) . Here we 
want to focus on general functional forms of the coupling strength 
in the context of coupled quintessence models where the interac- 
tion term is proportional to the energy density of the matter cou- 
pled fluid. One of the main motivations behind the introduction 
of a time dependence of the DE coupling to the matter sectors - 
besides the fact that an evolving coupling is per se a more gen- 
eral and natural assumption th an a constant in teraction strength 
- lies in the recent discovery JBaldi et alJuOlO.) that the effects 
of the DE-CDM interaction on the formation and evolution of 
structures at small scales, in particular in the nonlinear regime, 
might help alleviating the tensions between the ACDM model and 
a series of astrophysical observations. These range from e.g. the 



2003; Allen et al."2004l; [vikhlinin et al.l2006i ; lLaRooue et alj26oi 
.McCarthv et al.. .20071) . to the so called "cusp-core" probl em for 
the density profiles of the CDM halos of dwarf galaxies JMoord 
Ll994';'Flores & Primack"l994';'Simon et al."2003'), of spiral galax; 
ies (Navarro & Steinmetz 2000; Salucci & Burkert 2000; Salucci] 
2001; Bin nev & Evan s 2001), and of galaxy clusters JSand et 






.2002. . 20041; iNewmanetal ...2009) . or to the so called "Dark Flow" 



problem JWatkins et alj 120081) . Furthermore, some new potential 
challenges to the ACDM scenario have been recently reported 
base d on the detection of very massiv e high-redshift clusters (see 
e.g. |jee et al.ll2009l ; iRosati et alj|2009l) or on the observed dynam 



ical p roperties of CDM halo satellites JLee & Komatsull2010l ; lLei 

bold) . 



In their recent studv lBaldi et alj ( 120100 showed by means of a 
series of high-resolution N-body simulations how the DE-CDM in- 
teraction - for the case of constant coupling - could alleviate some 
of these problems; in particular, it was shown how the interaction 
can reduce the "cuspyness" of massive CDM halos, thereby go- 
ing in the right direction for a solution of the "cusp-core" problem. 
Nevertheless, the magnitude of this effect is strongly limited by the 
tight observational constraints on constant coupling models (as e.g. 
ISean et al.l2008l ; lLa Vacca et al.ll2009l) that put a firm bound to the 
maximum allowed value of the coupling. It is therefore natural to 
speculate about the possibility that a time dependent coupling with 
a large value during the late stages of structure formation and a 
progressively smaller value at high redshift might increase signifi- 
cantly the impact of the new physics introduced by the interaction 
on e.g. the density profiles of CDM halos without perturbing the 
overall evolution of the Universe beyond the present observational 
limits. 



The pre sent work constitu tes the natural extension of the anal- 
ysis done bv lBaldi et al.l ( 120101) to the case of time dependent cou- 
plings. In this paper we perform a complete numerical study of 
some quite general classes of coupling functions, starting from their 
background evolution up to the nonlinear regime of structure for- 
mation, and we present the first high-resolution hydrodynamical 
N-body simulations of structure formation in the context of inter- 
acting DE models with a time dependent coupling to date. 



The manuscript is organized as follows. In Section |2] we de- 
scribe the main features of the coupled DE models under investi- 
gation with a particular stress on the differences with respect to the 
standard constant coupling models that have been widely studied 
in the literature. More specifically, in Sec. l2.1l we discuss the back- 
ground equations and in Sec. |2.2| we illustrate the numerical meth- 
ods used to integrate such equations backwards in time. In Sec. 12.31 
we discuss observational constraints on the coupled DE scenario 
and a possible way to use the bounds derived for constant coupling 
models to check the viability of the variable coupling cosmologies 
investigated in the present work. In Sec. 12.41 we study linear per- 
turbations equations and we discuss the evolution of linear matter 
density fluctuations in variable coupling models. 
In Sec. [3] we briefly summarize the numerical methods used in the 
N-body simulations, and in Sec.|4]we present and discuss the re- 
sults of our runs. Finally, in Sec.|5]we draw our conclusions. 
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2 COUPLED DARK ENERGY COSMOLOGIES WITH 
TIME DEPENDENT COUPLINGS 

Coupled DE cosmologies have been widely studied in the last 
decade for what concerns their background and linear perturbations 
evo lution. In p a rticula r, we re fe r here to the derivations p resented 



lution. In p a rticula r, we re te r here to the derivations p r esente d 
Amendolal J2000l |2004) ; IPettorino & Baccigalupil ( l2008h : 



ISaldi et alj J201Ct) (BAIO, hereafter) and detailed references 
therein. Although these models of DE interactions have been 
mainly investigated for the simplified case of a constant coupling 
strength, the most natural scenario allows for a variation of the 
coupling along with the dynamic evolution of the scalar field 
(see e.g. lAmendolal l2004 h. and in the present work we want to 
investigate in detail this more general situation. To this end, we 
discuss here the most relevant features of coupled DE cosmologies, 
highlighting the main differences that arise for the case of time 
dependent couplings with respect to the standard case of a constant 
coupling. For this reason we will keep explicit the dependence of 
the coupling on the scalar field, and we will propose below some 
possible expressions for its functional form. 

We consider a DE model based on the dynamical evolution of 
a classical scalar field (ji rolling down a self-interaction potential 
Vi^tf))^ such that its intrinsic energy density and pressure can be 
expressed as: 

P<t> = 1:9^'' d^cj>ducj> + V((?!)) , 



P0 = \g^"d^^d,.^ - 1/(0) , 



(1) 

(2) 



where g^" is the metric tensor. The interaction of the scalar field 
with other fluids can be expressed as a source term in the conserva- 
tion equations for the different components of the Universe: 






VuT, 



(4'> 



-Q{t){4')T(i)Vv(t), 

^[Q(,)(0)T(,)]V.0, 



(3) 
(4) 



where Vp represents a covariant derivative, and T/^a^ is the stress- 
energy tensor ~ and Tf^j its trace ~ of the i-th component of the 
Universe, with i = c for CDM, b for baryons, 7 for radiation, n 
for neutrinos. Since the system of Eqs. ( I3l4t does not violate the 
conservation of the total stress-energy tensor of the Universe: 



'Z^'^(i)v 



0, 



(5) 



this type of interaction is consistent with general covariance and 
does not modify Einstein equations. 

It is also interesting to notice that radiation and relativistic 
neutrinos always remain uncoupled, since the stress-energy ten- 
sor of relativistic particles (subscript r) is traceless, T(^r) ~ 0. 
For massive neutrinos, instead, if Qn{<f)) 7^ the coupling term 
would become effective as soon as they become nonrelativistic, 
and such mechanism could provide solu tions to th e cosmic "co- 
incidence problem", as proposed in Amendola et alj (i2008}. 

A coupling of DE to baryonic particles is tightly c onstrained 
by S olar System tests of scalar-tensor theories (see e.g. lEllis et al.l 
[l989 : Carroll 1998; Will 2001). However, following the idea first 
proposed by Damour et al. ( 1990), one could consider models with 
different coupling strengths to baryons and CDM, for which the 
observational constraints become much broader. Therefore, since 
in this paper we are mainly interested in investigating the effects of 
a possible coupling between DE and CDM, we assume the baryonic 
coupling Q{b) {4>) to be vanishing at all times. 



To further simplify the system, we also assume neutrinos to 
be massless, such that they remain effectively uncoupled at all 
times and can therefore be considered part of the fraction of rel- 
ativistic particles of the U niverse, together with p hotons, such that 
Pr = p-y + pn- However, iLa Vacca et alj ( 120091) showed how re- 
leasing this assumption might have a relevant impact on the obser- 
vationally allowed range of coupling values, and we will discuss 
this possibility in Sec. 12. 31 

With all these assumptions we are left with only one non- 
vanishing coupling function, the DE-CDM coupling Qc[(l>), and 
the system of Eqs. OI4t . in a flat FLRW metric described by the 
line element [j: 



ds = —dt +a (t) iSijdx^dx-' ] , 



(6) 



results in the following set of dynamic equations: 



Pc + SHpc = — 
Pb + SHpb — , 

Pr + 411 Pr = , 



!*<«€ 



(7) 



3-*^ = T72- (pr +Pc+Pb+ P4>) , 

where an overdot represents a derivative with respect to the cosmic 
time t, H = a/a, Mpi = I/VSttG is the reduced Planck mass, 
and where we have followed the notation of iAmendola 12000,) by 
definingj: 



/3c(0) 



MpiQc{4>) ■ 



(8) 



An immediate consequence of the coupling source terms in 
Eqs. {Tjl is that the density of coupled matter species (CDM in our 



case) does not scale like a 
equation: 

pc{a) = pc{ao)a~ 



, but follows an evolution given by the 



/273//3e{0)d0/A/j,i 



(9) 



where the exponential extra factor accounts for the direct exchange 
of energy between DE and CDM. If one assumes that the number 
density of matter particles is conserved (i.e. particels are neither 
created nor destroyed), the direct consequence of Eqn. ^ is that the 
mass of coupled matter particles has to change in time according to 
the evolution of the scalar field: 



(10) 



M,{a) = Af,(ao)e-V2''3//3e(«d0/A/p, 



In this study we will always assume an exponential form for 
the scalar field self-interaction potential: 



V{(P) = Ae-^ 



/2/3a,p/Mpi 



(11) 



where a > is a free parameter of the model. We normalize the 
scalar field evolution such that (jf>(io) ~ 0, which implies A = 
V{to). 

Exponential potentials for the self interaction of a scalar 
fiel d have been propo s ed in the context of inflation by 
e.g. lLucchin & Matarresd ( Il985al lbl) and for the case of late-time 



^ We assume the speed of light c to be unity. 

^ Please notice that the definition of the coupling /3c in the present work is 

\/3/2 times the one adopted in BAIO. 
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acceleration by e.g. lWettericI] h988l) : iFerreira & Joyed ( Il998h : 
Ame ndola (20QO|). In both cases they determine the existence of 
attractor solutions that are almost independent from the scalar field 
initial conditions. 



2.1 Background evolution 

Following an d gen eralizing the approach devised by 
ICopeland et alj (1998), and subsequently largely applied to 
the s t udy o f coupled quintessence models (see e.g. lAmendolal 
bOOd. l2004h . we can rewrite the system of Eqs. ([7]( as an au- 
tonomous dynamical system by introducing the dimensionless 
variables: 



^/vW^ _ \fp^ _ \fp^ 



Mpi^/& 



MpiH 



MpiH 



MpiH 



(12) 

where now a prime represents a derivative with respect to the e- 
folding time N = ln(a). With these new variables the dimension- 
less density parameters Qi = pi/{3H'^AIpi) of the different cos- 
mic components can be expressed as: 



fib 



■ V ,0,,. ~ r , ^4> 



2 , 2 
X +y . 



(13) 



Furthermore, in this work we will always assume the Universe to 



be exactly flat, such that f^c = 1 — x — y 
recast the system of Eqs. l[7} in the form: 

x' — — (3x — Sy + r — 3) + ay 



y = 



, We can then 



+/?e(0) ( 


1- 


2 

X - 


2 2 

-y -r 


2 
— V 


- f(3^^- 


'3y^ 


+ 


r^+3)~ 


axy 


- ^(3-^- 


3y' 


+ 


r'-l), 




- i(3-^- 


3y' 


+ 


r'), 





(14) 



H' = 



H 



(3a;^ - 3y^ + r^ -f 3) 



The dynamical evolution of the system ( I14t has been thoroughly 
studied in the past, and analytic solutions for the critical points of 
the system and for their stability have been found for the case of 
a constant coupling function /3c(<^) = Pc (Amendola 200 0). One 
of the most prominent features of these cosmologies consists in 
the existence of a matter dominated scaling regime - which has 
been called a "(j>- Matter Dominated Epoch" (^MDE, hereafter) in 
lAmendolal J200G) - between the DE scalar field and the coupled 
matter component, where the two fluids keep a constant ratio of en- 
ergy densities before the final accelerated attractor is reached. Here 
we generalize the system ( I14t to the case of non-constant couplings 
investigated in the present work. For variable couplings, in fact, we 
need to rewrite the coupling function I3c{(t>) in terms of the dimen- 
sionless variables l ll2t in order to close the system. 

In this work, we will consider three possible different forms 
for the time evolution of the coupling function /3c (0) which are 
discussed below, and for each of these forms we will solve numer- 
ically the system il4t to get the background evolution of each spe- 
cific model. 



2.1.1 Coupling proportional to a power of the scale factor 

We start exploring a rather phenomenological form for the time 
evolution of the coupling strength, where I3c{4>) evolves as a power 



of the cosmological scale factor a: 



(15) 



This is a completely phenomenological parametrization of the time 
evolution of the coupling, since it does not depend directly on the 
dynamics of the scalar field (j), and might therefore look quite un- 
physical, but is particularly easy to implement and integrate, and 
can already give an idea of the main features that characterize the 
dynamics of a variable coupling model. For this case, in order to 
close the system l ll4t . we can just substitute /3c (0) with its ex- 
pression in terms of the dimensionless time variable of our system, 
which is the e-folding time A'^, as: 



/3c(9i) ^ Poe 



/3lJV 



(16) 



2.1.2 Coupling proportional to Q,p 

Another phenomenological possibility for the time variation of the 
coupling is to relate the coupling strength to the fractional DE den- 
sity during cosmic evolution. In this case, the evolution of the cou- 
pling depends on the dynamics of the scalar field as: 

n<4 



Pc{cb) = l3o- 



(17) 



and one can implement this type of coupling in the system J14t just 
by substituting /3c{<j)) as: 






(18) 



where we have indicated with a subscript the value at the present 
time N = 0. 



2.1.3 Exponential coupling 

The most natural form for the evolution of the coupling is given 
by a direct dependence of /3c on the value of the scalar field (j), as 
we stressed above. In particular, we consider here an exponential 
coupling in the form: 



PcW 



(19) 



Exponential forms of the c oupling strength b etween DE and CDM 
have been proposed in e.g. jAmendolal (2004). For this choice, a di- 
rect implementation of the coupling function /3c (0) into the system 
l ll4t is not possible since none of our dimensionless variables l ll2t 
represents the scalar field itself. 

One possible way to include the coupling l |19l l into our system 
of dimensionless equations consists in expressing /3c (0) in terms 
of the potential variable y, such that: 



/3c(</>) 



H'y' 



/3o 



ioVo 






/3/2/3i/q 



(20) 



An alternative possibility is to add a further equation to the 
system ( I14t . and to treat the scalar field as an independent degree 
of freedom. If we define 



e = ^/Mpi , 
from the definition of the variable x we get that 

C' = V6x . 



(21) 



(22) 
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By including Eqn. M2\ into the system l ll4t we can then express 
the coupling function Pc{(l>) as: 



/3.(<^) = /Joe'' 



(23) 



and close the system. 

These two approaches are both possible, but we found the lat- 
ter be simpler to handle numerically, and therefore we will prefer it 
for the numerical integration of the cosmological background evo- 
lution equations described in the next section. 



2.2 Numerical integration of tlie bacl^ground equations 

In order to study the background evolution of cosmological mod- 
els with a time dependent coupling between DE and CDM, we in- 
tegrate numerically the system of ordinary differential equations 
l ll4t plus the additional equation J22t for a series of possible mod- 
els with the three different types of coupling evolution discussed in 
Sec. 12.11 Since we want to find viable solutions for the expansion 
history, i.e. we want to realize a cosmological background evolu- 
tion that ends up at the present time with the observed values of 
the cosmological parameters, we integrate the system backwards in 
time, assuming the fiducial values of the cosmological parameters 
at z = as our initial conditions. This procedure is not straightfor- 
ward, and deserves to be briefly discussed. The main problem relies 
in the fact that the stability of the critical points of the autonomous 
system l ll4t is inverted under time inversion: stable points become 
unstable and vice versa. This makes the system be attracted towards 
the "wrong" solution in the backwards integration. 

To overcome the problem, we have devised an algorithm that 
starting from the desired cosmological parameters at 2 = inte- 
grates the equations backwards in time, and at the end of the inte- 
gration checks whether the system has reached the desired solution, 
which is given by a radiation dominated phase (r ~ 1). If the cor- 
rect solution has not been found, the initial conditions at 2; = are 
adjusted according to some specific prescription and the system is 
integrated again. This procedure is repeated until the correct solu- 
tion is found. The adjustment of the initial conditions for the system 
at 2 = consists in changing the ratio of kinetic to potential energy 
for the scalar field (j), consequently changing the present value of its 
equation of state parameter w^, defined as: 



w^ 



P4, 



y 



x^ + y-^ 



(24) 



Therefore, all our final "correct solutions" will share the same cos- 
mological parameters at z = except for w^ which is allowed to 
vary in a range [—1.0, —0.9]. 

For our integrations we assume the most updated set of max- 
imum likelihood cosmological parameters of WMAP-7 combined 
with BAO measurements and the Hubble Space T elescope (HST) 
determ ination of the Hubble constant as reported in lKomatsu et al.l 
( l2010h . which are listed in Tabled 

In order to check the convergence of our method we first apply 
our integration algorithm to the fiducial ACDM model by imposing 
Wcf, = — 1 (i.e. xo = , yo — V^de), and then we extend it to the 
known cases of an uncoupled scalar field, and a few coupled scalar 
field models with constant couplings. The numerical solutions of 
the backwards integrations for these well known models reproduce 
all the predicted features of uncoupled and coupled quintessence 
models respectively, like e.g. the existence of a "<j!>-MDE" phase 
for the models with a constant coupling. Finally, we can use our 





Parameter 


Value 








Ho 70.2kms-iMpc- 


-1 






f^CDM 


0.227 








f^DE 


0.728 








CS 


0.807 








nt, 


0.0455 








ris 


0.961 






Table 1. Cosmological parameters for our set of models, consistent with the 


WMAP 7 year 


results for a ACDM cosmoloev tKomatsu et alj2010iy 


Model 


Type of coupling 


/3o 


/3i 


w^{z = ff) 


ACDM 


- 


- 


- 


-1.0 


EXPOOO 


- 


- 


- 


-0.999 


EXP002 


M<l>) = Po 


0.1 


- 


-0.995 


EXP005 


" 


0.25 


- 


-0.984 


EXPO 10 


" 


0.5 


- 


-0.945 


BEOS 


" 


0.067 


- 


-0.997 


LV09 


" 


0.17 


- 


-0.991 


EXPOlOa 


/3c(</.) = Poa^^ 


0.5 


1.0 


-0.981 


EXP010a2 


" 


0.5 


2.0 


-0.990 


EXPOlOaS 


" 


0.5 


3.0 


-0.993 


EXPOlSa 


" 


0.75 


1.0 


-0.963 


EXP015a2 


" 


0.75 


2.0 


-0.981 


EXP015a3 


" 


0.75 


3.0 


-0.988 


EXP020a 


" 


1.0 


1.0 


-0.938 


EXP020a2 


" 


1.0 


2.0 


-0.970 


EXP020a3 


" 


1.0 


2.0 


-0.981 


EXPOIODE 


/3c(0) = /3on<^/f20,o 


0.5 


- 


-0.986 


EXP015DE 


" 


0.75 


- 


-0.975 


EXP020DE 


" 


1.0 


- 


-0.960 


EXPOlOe 


/3c(</') = /3oe''i'^'"'^ 


0.5 


1.0 


-0.966 


EXP010e2 


" 


0.5 


2.0 


-0.972 


EXPOlOeS 


" 


0.5 


3.0 


-0.976 


EXPOlOe 10 


" 


0.5 


10.0 


-0.987 


EXPOlOe 15 


" 


0.5 


15.0 


-0.989 


EXPOlSe 


" 


0.75 


1.0 


-0.935 


EXP015e2 


" 


0.75 


2.0 


-0.952 


EXP015e3 


" 


0.75 


3.0 


-0.960 


EXPOlSe 10 


" 


0.75 


10.0 


-0.979 


EXP015el5 


" 


0.75 


15.0 


-0.984 


EXP020e2 


" 


1.0 


2.0 


-0.926 


EXP020e3 


" 


1.0 


3.0 


-0.941 


EXP020el0 


" 


1.0 


10.0 


-0.971 


EXP020el5 


" 


1.0 


15.0 


-0.976 



Table 2. List of all the cosmological models considered in the present back- 
ground evolution analysis. The models are divided based on the type of cou- 
pling function fic{4') between the DE scalar field and the CDM fluid. The 
last column indicates the value of the equation of state of the scalar field 
ui0 at 2 = obtained from the numerical integration of the background 
evolution, and represents the closest value of w^ to the cosmological con- 
stant value of — 1 that is possible to obtain for each model and for the set of 
cosmological parameters listed in Table[T] The potential slope a is equal to 
0.1 in all the cosmologies. 



algorithm to integrate a number of variable coupling models, and 
obtain for each of them the detailed background evolution. 

All the models, with the specific values of the coupling pa- 
rameters /?o and /3i , and the value of the resulting equation of state 
parameter w^ at 2 = are listed in Table [2] The potential slope 
parameter has the same value a = 0.1 in all the cosmologies. 
The evolution with redshift of the mass correction factor given by 
Eqn.[TO]is shown in Fig.[T]for all the models listed in Table|2] 

Even in the absence of analytic solutions for the cosmologi- 
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Particle mass evolution for Po = 0.5 




Particle mass evolution for P(, = 0.75 




Particle mass evolution for Po = 1 .0 




Figure 1. Mass correction of CDM particles as a function of redshift for 
the different interacting DE models under investigation. The three panels 
refer to models with three different values of the coupling at z = 0, respec- 
tively f3c{4>o) = 0.5 (upper panel), /3c(0o) = 0.75 (middle panel), and 
Pc{<j>o) = 1-0 (lower panel). The color legend for all the models is given 
in the upper panel 



cal background evolution of variable coupling models, our nume- 
rical integrations allow to identify some relevant features of these 
cosmologies and to compare them with the familiar case of con- 
stant couplings. In particular, it is very interesting to notice how 
the time variation of the coupling affects the background evolution 
during matter domination. As already mentioned above, a constant 
coupling between DE and CDM gives rise to a metastable scaling 



solution during matter domination where the scalar field (j) and the 
coupled matter fluids (CDM in our case) evolve with a constant ra- 
tio of energy densities, called "<j!>MDE", which represents one of 
the most peculiar features of the coupled DE scenario, besides pro- 
viding a way to put constraints o n the models due to the presence of 
a sizeable fraction of early DE (Wette richll2004l : lDoran et al.ll200ll : 
[Doran & Robbers 2006) at high redshifts. 

For the idealized case of a Universe containing only DE and 
CDM interacting with a constant coupling /3c, the DE fractional en- 
ergy density Q.^ during the "(^iMDE" phase has the constant value: 



n4(/.MDE) = ^/3? . 



(25) 



In our variable coupling models, where the coupling is large at 
the present time and progressively decreases at higher redshifts, 
we find a quite different picture for the background evolution. In 
particular, there is no "<j!>MDE" solution in our cosmologies, and 
no scaling regime between the scalar field (j) and CDM is found. 
Among the different models investigated, we can distinguish two 
quite different behaviors for the different types of couplings consid- 
ered in our integrations. The situation is illustrated in Fig.[2l where 
the time evolution of the scalar field fractional energy density Q.^ 
and the redshift evolution of the scalar field equation of state pa- 
rameter w^ are plotted for some of the models under investigation. 
For an easier readability of the plots, and with no loss of gener- 
ality, we show in Fig. [2] only the models that will be used for our 
N-body simulations since they clearly show the different impact on 
the background evolution between the exponential coupling mod- 
els (/3c (0) oc e'^i'*/") and the scale factor models (/3c (0) oc a^^). 
The remaining class of couplings 0c{(l>) ex: fi^) is found to have 
a very similar behavior to the latter one. In addition, we also plot 
for comparison the same quantities for one of the constant coupling 
scenarios studied in BAIO, the RP5 model, which features a cou- 
pling of ^Sc = 0.25. 

In the left panel of Fig. |2] the solid lines represent the evolu- 
tion of the DE fractional energy density fl^. The "(/iMDE" phase 
is clearly visible for the constant coupling model RP5, while none 
of the variable coupling models follows a similar evolution. How- 
ever, while the EXP010a2 and EXP015a3 models, where the cou- 
pling scales like a power of a, present a constant decay with red- 
shift of the DE fractional energy density, with no significant dif- 
ference from a ACDM evolution, the exponential coupling mod- 
els EXP010e2, EXPOlOeS, and EXP015e3 show an intermediate 
behavior, which we call a "Growing </)MDE" solution, where the 
fraction of DE ^^, although not constant during all matter dom- 
ination, evolves significantly slower than in the ACDM scenario. 
In order to compare with the constant coupling case, we have also 
plotted (as dotted lines) the analytic solution for f20(0MDE) given 
by Eqn. l |25t . derived for the case of constant /3c, also for the vari- 
able coupling models. The gap between the theoretical value and 
the actual evolution of f2</, for the constant coupling RP5 model is 
due to the inclusion in our integrations of the uncoupled baryonic 
component that perturbs the solution given by Eqn. l |25l l for the 
idealized case of Qt ~ 0. 

It is very interesting to notice that during the "Growing 
(;/)MDE" phase, the models with exponential couplings follow 
the same evolution given by the extension of the solution l |25l l 
to the case of a variable coupling /Jc(0), with a comparable gap 
as for the RP5 model due to the presence of uncoupled baryons. 
The denomination "Growing (jiMDE" therefore seems appropriate 
since during this phase the system evolves as for the case of the 
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Figure 2. Left Panel: The evolution of the DE fractional energy density Q,^ (soHd lines) for a few models of variable coupHng, and in addition for the constant 
coupling models RP5 studied in BAIO. The "</<MDE" phase is clearly visible for the RP5 model, and is completely absent for the phenomenological coupling 
models EXP010a2 and EXP015a3, that show no significant difference with respect to a ACDM evolution. The the exponential coupling models have the 
intermediate behavior of a "Growing (/)MDE". Dotted lines represent the value of the analytic "t^MDE" solution I Amendola 2000) extended to the case of 
growing couplings. Right Panel: Equation of state parameter w^ for the same models as in the left panel. The phenomenological coupling models have always 
an equation of state parameter very close to the cosmological constant case, while the exponential coupling models follow a similar behavior to the constant 
coupling model RP5. 



standard "(/)MDE" but with a growing value of the coupling. 

The reason for the strong quahtative difference in the back- 
ground solutions found for the different classes of coupling func- 
tions can be understood by looking at the right panel of Fig. |2] 
where the DE equation of state parameter Wtt, is plotted as a 
function of redshift. While the scale factor models (EXP010a2, 
EXP015a3) always have an equation of state parameter very close 
to —1, the exponential coupling models show a significant variation 
of w^ with redshift, with a similar trend to the constant coupling 
model RP5. 

This strikingly different behavior of the equation of state 
clearly shows that the former class of models requires the same 
level of fine tuning of the ACDM concordance cosmology, and 
therefore cannot be expected to provide any dynamical solution to 
the DE "fine tuning problem". On the contrary, the "fine tuning 
problem" might still be alleviated by the latter class of models, due 
to their sensibly slower evolution of the scalar field energy density 
with respect to ACDM, thereby retaining one of the main motiva- 
tions that are behind the whole interacting DE scenario. 

However, as we will discuss in detail in Sec. [4] the different 
dynamics of these two types of variable coupling models is not 
relevant only for the cosmic background evolution, but will also 
have a strong impact on the structure formation processes taking 
place in the context of these different cosmologies, in particular for 
what concerns the highly nonlinear regime of structure formation. 



2.3 Observational constraints and model selection criteria 

The observational constraints on the DE-CDM interaction have be- 
come ever tighter in the last decade. This is mainly due to the 
inclusion of low redshift probes and Large Scale Structure (LSS) 
data in the analysis rather than to an improved quality of CMB 
data, which have been for a long time the main, if not the only, 
source of constraints for interacting DE models. In fact, the most 



stringent constraints to date on the coupling strength for constant 
coupling models have been presented in Bean et al. (2008) (BEOS, 
hereafter) and come from a combined analysis of high and low red- 
shift probes of the expansion history of the Universe as CMB data, 
HST measurements of the Hubble constant Hq, Baryon Acoustic 
Oscillations data, Snia measurements of the low redshift expan- 
sion history, and Sloan Digital Sky Survey (SDSS) measurements 
of the matter power spectrum at low redshifts. By combining all 
these datasets, BEOS came up with a constraint of |/3cj < 0.067 
at the 95% C.L. for a constant coupling interacting DE scenario. 
However, it is interesting to notice that the same authors quote a 
limit of |/3ci < 0.13 from the analysis of CMB data alone, which 
is no much tighter constraint than the upper bound of |/3cl < .15, 
still based on CMB data only, reported bv lAmendola 1 20001) al - 
most a decade before . This shows, as confirmed by e.g. lXial ( 120090 : 
IValiviita et al.l ( 120101) , that CMB data alone are not able to constrain 
much further the value of the coupling, and the inclusion of low red- 
shift probes as well as LSS data is essential to tighten the bounds 
down to the accuracy reported by BEOS. In addition to these com- 
bined constraints a new and completely independent bound on the 
coupling based on the comparison of detailed hydrodynamical N- 
body simulations with the observed propertie s of Lyman-ct absorp- 
tion systems has been recently presented by iBaldi & Viell ( 120100 . 
with a 2-0 limit on the coupling of /? < 0.15 based on Lyman- 
Q observables only. 

An interesting discovery has nevertheless been recently re- 
ported bv lLa Vacca et alj ( l2009l) (LV09, hereafter), that found a de- 
generacy between the coupling amplitude and the average neutrino 
mass Ml,, showing how a value of M^ ~ 1 eV could broaden the 
allowed range for the coupling up to | /3c | < 0.17. Therefore, by 
dropping the assumption of massless neutrinos that we adopted in 
Sec.|2l it is possible to allow larger coupling values without running 
into conflict with present observational constraints on the back- 
ground expansion history. The inclusion of massive neutrinos with 
an average mass up to M^ ~ 1 eV would not significantly change 
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the results we have derived so far for the background evolution as 
long as the neutrinos remain uncoupled, i.e. I3n{<j}) = 0. 

We will assume here the values of /3c = 0.067 reported by 
BEOS and of /3c ~ 0.17 derived by LV09 as two limiting values for 
the allowed coupling in a constant interaction model, for the cases 
of massless and massive (M,j ~ 1.0 eV) neutrinos, respectively. 
These values correspond to the BEOS and the LV09 models in Table 

HI 

Some attempts to put observational constraints on more gen- 
eral classes of interacting DE models, inclu ding some form s 
of va riable couplin gs, have been present ed in ICai & Sj ( 120100 : 
iGuoetal. (200 7); C osta & Alcanizl ( l2010h and further generalized 
bv IWeil ( 120101) . However, such constraints are not directly appli- 
cable to the scalar field models under investigation in the present 
work, due to the assumption of a constant DE equation of state pa- 
rameter that is made in the analysis of all the cited authors. Also, 
the fraction of uncoupled baryonic mass should always be included 
in the analysis in order to derive meaningful constraints. 

Therefore, a full likelihood analysis of variable coupling mod- 
els against presently available data is still missing, and will be 
clearly necessary to fully assess the viability of these models and 
put constraints on the coupling function. However, such an effort 
goes beyond the scope of this paper, and we leave it for future work. 
Here we just assume a simple criterion to test and select our mod- 
els, based on a direct comparison of their background evolution 
with the one derived for our assumed limiting models BEOS and 
LV09. We compute for all our models the evolution of the luminos- 
ity distance: 



dL{z) 

between z = and z 
diameter distance: 

dA{z) = 



l + z 
Ho 



Hodz 
lliz) 



(26) 



3, and the evolution of the angular- 



Hpdz' 

{l + z)HoJo H{z) 



(27) 



up to last scattering surface (z ~ 1100), and we compare their 
relative evolution with respect to the fiducial ACDM model to the 
limiting cases BEOS and LV09. The change in the angular-diameter 
distance at high redshift would affect CMB measurements, while a 
change in the evolution of the luminosity distance with redshift at 
small z would be constrained by low redshift data as e.g. Snia. 

Therefore, since most of the present bounds on interacting DE 
models are based on observations of the background expansion his- 
tory, we assume as a tentative selection criterion that any model that 
lies in between the fiducial ACDM cosmology and one of the two 
limiting scenarios BEOS and LV09 for what concerns both its lumi- 
nosity distance at low redshifts and its angular-diameter distance at 
high redshifts should be considered compatible with observations 
and accepted as viable, while it should be rejected otherwise. The 
situation is illustrated in the two plots shown in Fig. [S] where the 
ratio of luminosity distance (left panel) and angular-diameter dis- 
tance (right panel) for all our models with respect to ACDM is plot- 
ted against redshift. The two limiting scenarios we are considering, 
BEOS and LV09, are shown as the two thick solid lines, and the 
allowed regions for the two limits are represented by the dark-grey 
and light-grey shaded areas, respectively. 

As Fig. [3] shows, while only a few of our models pass the 
test of the most stringent constraint given by BEOS, most of them 
are found to be compatible with the other limiting case LV09. It 
is nevertheless interesting to notice that none of the exponential 
coupling models seems to be consistent with the limit given by the 



analysis of BEOS, for which they would probably require a much 
larger value of the exponential coupling slope /3i . 

The selection criterion presented here is clearly not rigorous 
enough to be taken as a definitive check of the viability of the 
models under investigation, but we find it a simple and useful 
guideline to select on which models we should invest our com- 
putational time for the high-resolution N-body simulations de- 
scribed in Sec. [3] Based on this criterion we therefore decide to 
run high-resolution simulations for the two models consistent with 
BEOS, namely the EXP010a2 and the EXP015a3 models, and for 
three of the more physically motivated exponential coupling mod- 
els that passed at least the test of the bounds of LV09, namely the 
EXP010e2, EXPOlOeS and EXP015e3 models. 

The choice of the latter models among all the models com- 
patible with LV09 is not random. In fact, although all the models 
considered in the present work are specifically built in order to have 
very large values of the coupling in a redshift range relevant for the 
nonlinear stages of structure formation, and are therefore expected 
to imprint sizeable features in the properties of highly nonlinear 
structures, the way in which the coupling affects the newtonian dy- 
namics of CDM particles does not depend on the coupling strength 
alone, but also on the evolution of the DE internal degrees of free- 
dom like the scalar field velocity cj), as we will discuss in detail 
in Sec. 12.41 The models we have chosen for our N-body treatment 
feature substantially different behaviors of the scalar field veloc- 
ity (j> (the DE kinetic degree of freedom x in our dimensionless 
framework) which will play a significant role in determining how 
the DE-CDM interaction can influence the properties of collapsed 
structures like e.g. the halo density profiles, as we will discuss in 
full detail in Sec. [441 



2.4 Linear perturbations 

We now study the evolution of density perturbations according to 
linear perturbation theory for the variable coupling models intro- 
duced above. This is an interesting issue on its own since such anal- 
ysis has not been done previously for the case of variable couplings. 
Linear perturbations equations in the context of coupled DE mod- 
els have been derived in the literature (see e.g. Am endolal (2000|, 
2004); Pettorino & Baccigalupi (200S), BAIO), also for the case of 
a scalar field dependent coupling, but to our knowledge the solution 
of such equations for variable coupling models has not been pre- 
sented before. We refer the interested reader to the cited literature, 
and references therein, for a detailed derivation of the perturbed 
equations, and we limit our discussion here to the main features of 
the perturbations growth that are relevant for the models consid- 
ered in this work. To ease the readablity of lengthy equations, in 
this section we do not keep explicit the dependence of the coupling 
on the scalar field. 

The direct solution of the linear density fluctuations evolution 
is also necessary in order to set up the initial conditions for our 
set of N-body simulations, since we want to be able to normalize 
our cosmological models to the same parameters at 2 = 0, which 
means that we also want all the cosmologies to end up with the 
same as, and therefore initial conditions have to be properly scaled 
with respect to each other according to the different values of their 
growth factor at the starting redshift of the simulations, which we 
assume to be 2 = 60 for all of our runs. 

We consider the perturbed metric in the longitudinal gauge 
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Figure 3. The ratio of the himinosity distance between 2 = and z = 3 {left panel) and of the angular-diameter distance between 2 = 10 and z = 1000 (right 
panel) to the respective ACDM evolution as a function of red shift . The two thick solid b lack lines in the two plots represent the two limiting cases considered 
in this work, i.e. the constraints given bv mean et alj (20081) and iLa Vacca et al J 420091) . The dark-grey and light-grey shaded regions that lie, respectively, 
between the ACDM value of 1 and each of the two limiting scenarios represent the allow ed regions accordin g to the selection criterion described in the text. 
As it can be noticed, only a few models are found to be compatible with the constraints of lBean et alj 420081) . while a large fraction of the models considered 
in the present work are compatible with the bounds of lLa Vacca et alj 42009h . 



and in the absence of anisotropic stress, given by: 



ds 



\~{l + 2^)dT'^ + {l-2^)5ijdx'dx^'\ , (28) 



where r is the conformal time and <I> is the gravitational potential. 
The conformal Hubble function is given by H = {da/dT)/a — 
aH. Then, one can define the following perturbation variables: 



5c = 5pc/ pc , 5b = 5pb/pb : 

Ui = adxi / {T-Ldt) , V ■ Uc = 
if = 5(t)/{Mp,V6). 



9c,V -Ub- 



(29) 
(30) 
(31) 



In Fourier space we can also define the scale parameter A = T-L/k. 
Since we are interested here in the evolution of matter density per- 
turbations during matter domination, we discard for simplicity the 

radiation density fluctuations. 

Finally , we define the dimensionless scalar mass as in I Amendolal 
( |2004|) : 



ml 



1 dV((^) 



H H 



o 2 2 

2a y . 



(32) 



With all these definitions, the perturbed evolution eq uations for 
each component in Fourie r space can be written as (see I Amendolal 
l2004l : |Pettorino & Baccig alupi 2008, for a complete derivation): 



5' 



5'b 



1 + 



Oc + 3$' - 2l3cip' - 213'cip , 

— -^fl„_r 1 ^„ 4- \~^' 

H 
+ 3$', 



2l3cx]ec + X~''{^-2l3c<fi)., 



1 + 



H 



+ >-"$, 



+ 



n 



v' + (a- 



-, 2 



«c/3c 



-4$'s - 2y^a$ = /3cf^e(5c + 2-1') , 
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(33) 
(34) 
(35) 
(36) 

(37) 



As it has been noticed in I Amendolal d2004l) . it is very inter- 
esting to point out that the variation of the coupling contributes to 
the equation for the evolution of the scalar field fluctuations ( 1371 ) 
as an effective mass term. However, we want to focus our attention 
here on a point that has not been sufficiently stressed in previous 
literature, namely the fact that this effective "coupling mass" 



.- 2 



^AW 



(38) 



can take both positive and negative signs depending on the signs 
of X and /?c(0)- As a particular case, if one assumes a direct de- 
pendence of the coupling on the scalar field, as it is the case for 
the exponential coupling introduced in Sec. 12. 1.3] the sign of m|^ 
is fully determined by the derivative of the coupling function with 
respect to the scalar field itself: 



m^^ — v&Q,, 



dPcW 



(39) 



It is therefore evident that a negative effective mass rh = rh^ — 
m^^ could arise for couplings that grow in time, as it is the case for 
the models considered in the present work. If this happens, large- 
scale instabilities might arise in the growth of density perturbations. 

For the moment we ignore this potential instability issue by as- 
suming that in the newtonian limit (small scales, A ^ 1) that is rel- 
evant for our N-body simulations, the total mass term in Eqn. l l37t 
can be discarded as compared to the A~'^ term. This is true as long 
as \fr?\ ~ C(l), which will be tested later on for all the specific 
models under investigation. 

In any case, the issue of potential large-scale instabilities due 
to the coupling growth should be carefully investigated in order to 
ensure the viability of the models, and w e leave such analysis for 
future work ( Amend ola & Baldillin prep J) . 

By applying the newtonian limit A ^ 1, the perturbations 
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equations for baryons and CDM become: 



(40) 



1 + ^ - 2/?cS 



Oc-^ l^bSi + nj,r,] , (41) 



SL 



i + ^^Ob-^lfttSt + ncSc] 



where we have defined: 



r, = 1 + 



4 I3!{^) 
3 1 + A2m2 



(42) 
(43) 

(44) 



Since from the newtonian hmit of Eqn. J37t one gets tha t 
the s calar field fluctuation ip y X (see again I Amendolal ( 120040 : 
IPettorino & Baccigalupi (2 0081) for a complete derivation), the term 
in P'^ in Eqn. ( I40t can be dropped as long as l3'^{(j)) ~ 0(1). Sim- 
ilarly, the term A m at the denominator in the second term of 



0(1), in which case 



(45) 



Eqn. ( I44t can be discarded as long as \m~\ 
one gets: 

Tc = 1 + ^p!(<p) . 

As we will show in Sec. 12.51 the former condition is always ful- 
filled for all the models under consideration, and we therefore as- 
sume it to hold. The latter condition, instead, due to the factor 1/x 
in the definition of the coupling mass m|^ (Eqn. I38t , might not 
be verified at all subhorizon scales for all our models, and has to 
be carefully evaluated case by case. In particular, as we will show 
later on, some of the models under investigation feature a slow-roll 
regime of the scalar field during matter domination, which makes 
the factor 1/x significantly large, with a corresponing increase of 

Assuming the validity of the condition P'^j rj}) ~ 0(1). the 
perturbations equations take the well known form d Amendolal 200CI . 
[2004; Pettorino & Baccig alupi.i2008,): 



Si 



SL 



1 + ^ - 2l3,x 



(46) 

ec-^ [njt + iicS^r,] , (47) 



i + j^]eb-^[ntSt + nj,] , 



(48) 
(49) 



which leads, by derivation of Eqs. l l46l48t , to the dynamic equa- 
tions for density fluctuations: 

<5" + (l + ^ - 2/3cS j S',-^ [ntSt + ftj.r,] = ,(50) 



5'^ +{i + ^)si-^[ntSt + ncSc] = o, 



(51) 



and to the vectorial acceleration equations in real space for baryon 
and CDM particles, that was first derived in its full generality in 
BAIO: 



Vc ~ —HVc 



Vb 



-Hvc 



E 



E 



GAf,(<^) ^-^ GMj 



Ti 



E- 



(52) 



(53) 



In Eqs. ( |52|53l l H = ii'[l+2/3c(<?!))a;] and d; is the peculiar velocity 
of a baryonic (subscript b) or a CDM (subscript c) test particle. 



Vi = axi. To obtain Eqs. ( 152153b we have discretized the mass 
distribution in the Universe by assuming: 



^cSc ~ y 

i^bSb = y 

j=b 



3Wa 



(54) 
(55) 



where 5{ri) stands for the Dirac distribution, Vi is the position of 
each matter particle with respect to our test particle, and the sum 
has to be intended as running over all the other matter particles in 
the Universe. 

Eqs. ( I50l51t are the equations that we have to solve in order 
to compute the specific growth factors that we need for setting up 
the initial conditions of our N-body simulations, while Eqs. ( |52|53l l 
are the modified newtonian dynamic equations that have to be im- 
plemented in the N-body algorithm in order to follow correctly the 
dynamics of particles in a coupled DE cosmology. 

It is importan t to notice at this stage, as it was shown in 
iMaccio et alj ( 120041) and BAIO, that the acceleration equation for a 
CDM particle is modified in three ways, due to the presence of a 
DE-CDM coupling. 

First, the friction term Hvc is modified by the presence of an 
extra contribution 2l3c{4>)xVc- It is interesting to notice, as antic- 
ipated in Sec. 12.31 that this extra friction does not depend on the 
coupling strength alone, but is modulated by the scalar field kinetic 
energy x. Therefore, the impact of the friction term on the dynam- 
ics of CDM particles will depend substantially on the background 
evolution of the scalar field, which will have a strong influence on 
how the coupling affects the properties of nonlinear structures. In 
Fig. |4] the coupling function /3c (</>), the kinetic energy x, and the 
friction term coefficient 2/3c((/))a; are plotted as a function of red- 
shift (panels a, b, and c, respectively) for all the models we have 
selected for simulations and additionally for the RP5 model stud- 
ied in BAIO, which has a constant coupling of /3c — 0.25. It is very 
interesting to notice how, despite the large values of I3c{(t)) for all 
the variable coupling models at low redshift, the magnitude of the 
friction term can be strongly suppressed in some of the models by 
the faster decay or the lower initial value of the DE kinetic energy 
X. This means that the friction term becomes more efficient the 
more the DE equation of state w^ departs from the cosmological 
constant value of —1. 

Second, the mass of CDM particles that generate the gravita- 
tional potential in which every particle is moving changes in time, 
due to the evolution of the scalar filed, according to Eqn. dlOt . This 
variation of CDM particles mass affects also the dynamics of the 
uncoupled baryonic particles, as it can be seen in Eqn. l |53| l, which 
will feel a decaying gravitational potential due to the CDM mass 
loss. 

Finally, the gravitational acceleration of CDM particles in- 
cludes an extra factor Pc accounting for the fifth-force mediated 
by the scalar field. So far, we have left this extra factor in its full 
generality, but for the set of N-body simulations carried out in this 
work we will limit the analysis to the case where the effective 
scalar mass m? is negligible, such that Vc is given by Eqn. ( 145 1 . 
This corresponds to a long-range scalar fifth-force, and we will dis- 
cuss in Sec. 12.51 under which circumstances this approximation is 
vali d. It is neverthele ss interesting to notice that, as discussed in 
e.g. Amendolal ( l2004r) . the more general case where the scalar mass 
is not negligible, given by Eqn. ( I44t . determines a scale dependence 
of the fifth-force interaction between CDM particles corresponding 
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Figure 4. This figure shows the evolution as a function of redshift of the coupling strength Pc{<t>) {panel a), of the scalar field velocity x as defined in Eqn. ll2l 
{panel b), and of the friction term coefficient given by 2a;/3c(</>) {panel c), for the five coupled DE models studied with our high-resolution N-body simulations 
and in addition for the constant coupling model RP5 studied in BAIO which is overplotted for comparison. It is worth noticing how the dynamic evolution of 
the scalar field modulates the effect of the coupling in the friction term: despite a large value of the couphng /3c (</>) at low redshifts the friction term is strongly 
suppressed in the models which feature a too fast decrease of the coupling with redshift due to the absence of a proper "</<MDE" or of a "Growing 0MDE" 
phase. 



to a Yukawa potential: 

4(r) 



GM 



l + l/3c^(0)e-*«'- 



Coupling derivatives 



(56) 



Cosmological models that feature only this type of modified 
gravitational potential in the dark sector have been studied in the lit- 
erature assuming differen t values of the characteristic mass scale m 
JGubser & Peeblesll2004l : JFarrar & Peeblesl |2004 JFarrar & Roser] 
120070 . a nd tested with numerical N-body simulations both at 
galactic JKesden & Kamio nkowski 2006; Keselman et al. 2 0091) 
and cosmological scales aNusser et al.l2005l : lHellwing et aU20100 . 
Simulations of massive scalar fields coupled to CDM, which fea- 
ture not only the Yukawa-like fifth-force of Eq. l |56| l, but also the 
other modifications of newtonian dynamics arising from the DE- 
CDM interaction which are described above, will be carried out in 
an upcoming publication. 



2.5 Scalar masses and stability conditions 

In the previous section we have encountered a series of conditions 
on the time variation of the scalar field coupling /3c (<?!>) that we as- 
sumed to be fulfilled in order to derive the final forms of the pertur- 
bations evolution equations l |50|51| l and of the particle acceleration 
equations I l52l53t within our coupled DE models. In this section 
we want to discuss in more detail such conditions and check their 
validity for the models under investigation. 

First of all, in deriving Eqn. l l46t we have assumed that 
Pc{4') ~ C'(l)- We will show here that this condition holds for all 
the models considered in this work. In particular, for the cases of 
an exponential coupling and of a coupling proportional to a power 
of the scale factor described in Sec. 12.1.31 and 12. l.Tl respectivelv. it 
is possible to estimate analytically the magnitude of /3c (0): 

/3c(</i) = Poa^' -^ P'M) = /3o/3ia''^ < /3o/3i , 
/3c(0) = /3oe'^i^/"^/3:(0) = ^a;/3o/3ie'5i^/". 



(57) 
(58) 



For the former case, one can easily compute that the largest possi- 
ble value for the models under investigation is given by /3c (0) ~ 3. 
For the latter, one can use the observation that in all the cosmolo- 
gies X < 0.15 to put an upper limit on the coupling derivative 
which is I3'^{4>) ~ 0.9 for the largest assumed values of /3o and /3i. 




Figure 5. The evolution of the coupling derivative ^'^{(p) as a function of 
the e-folding time N = Ina. The different colors represent the different 
types of time evolution of the coupling function f^c{<t>) while the different 
linestyles correspond to different values of the coupling at the present time. 
It is important to notice that none of the models presents a value of the 
coupling derivative larger than 0(1) at any stage of cosmic history. 



For the remaining case of a coupling proportional to the DE frac- 
tional density Q^ such analytic estimation is not possible, and we 
compute the derivative numerically. In any case, as shown in Fig.[5] 
the coupling derivative I3'^{(l>) is never larger than 0(1). Therefore, 
the simplification adopted in Sec. l2.4l is fully justified. 

Secondly, as we have seen in Eqn. l l37t . the variation of the 
coupling appears in the scalar field perturbations equation in the 
form of a "coupling mass" term mg2 , given by Eqn. i38t . such that 
the total effective mass takes the form: 



^ 2 



(59) 



This might lead to large-scale instabilities in the scalar field per- 
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turbations if the total effective mass is negative, wliich is the case 
for all the models under investigation in this work, as we will show 
below. Nevertheless, if the absolute value of the total effective mass 
|m^| is not too large, such instabilities might be confined to very 
large scales and the consequent growth of the large-scale gravita- 
tional potential might still be consistent with present observations. 
The total effective scalar mass for the three different types of cou- 
pling discussed in Sec. 12. li as a function of the e-folding time A'^ is 
shown in the three panels of Fig.|6l As the plots show, all our mod- 
els present a negative effective mass during all matter domination. 
Even though the absolute magnitude of the effective mass \rh?\ is 



not very large, a value of rh 



-10, that is realized for the more 



phenomenological models of /3c oc a^'^ and /3^ oc Q,^, might be- 
come relevant at scales of a few hundred Mpc. On the other hand, 
the more physically motivated exponential coupling models, as one 
can see in the right panel of Fig. |6l feature a much smaller am- 
plitude of the total effective scalar mass, whose absolute value is 
always \rh?\ ~ 0{1), and should therefore not show instabilities at 
any subhorizon scale. 

The issue of potential large-scale instabilities related to the 
behavior of the effective scalar mass in variab le coupling models 
will be d iscussed in detail in a dedicated work aAmendola & Baldil 
lin prepj) , but for the phenomenological study of variable coupling 
models addressed in this paper we assume the small scales of inter- 
est for N-body simulations (~ few tens of Mpc) not to be influenced 
by such potential large scale instabilities. Also, we will discard the 
Yukawa-like correction to the force law discussed in Sec |2.4| due to 
the small size of our simulation box, and we will always consider 
the scalar fifth-force to have infinite range. 



3 THE SIMULATIONS 

3.1 The simulations set 

The main focus of the present work is on the effects that coupled 
DE models with a time dependent coupling have on structure for- 
mation, and in particular on the properties of highly nonlinear struc- 
tures as cluster-sized or galaxy-sized CDM halos. This is done by 
means of a series of high-resolution hydrodynamical N-body simu- 
lations carried ou t with a suitably modified version of the numerical 
code GADGET-2 JSpringefcOOSh . 

N-body simulations of coupled DE models have been per- 
formed before for the case of a constant coupling by Maccio et alJ 
( 120041) and BAIO, and this work constitutes the natural extension 
to a time dependent coupling of the latter publication. In particular, 
we use here the same modified version of GADGET-2 that was used 
in BAIO, since the implementation presented and used in that work 
includes also variable coupling models, even though it had been 
used so far only for the simplified case of a constant coupling. 

To our knowledge, this work presents the first N-body simula- 
tions of coupled DE models with time-dependent couplings carried 
out to date. Previous attempts to use N-body simulations for other 
types of modified newtonian gravity (i.e. dif ferent from the cou - 
pled DE scenario) have been discussed e.g. in Nusser et alj 1 20050: 
IStabenau & Jain' ('200^;'Springel & Farrai' ('2007^;'Laszlo & Bear^ 
(2008); Sutter & Ricker (2008); Ovaizu (2008); Li & Zhao (2009); 
fHellwing et alj ( l2010l) ; lDe Boni et alj ( |201 0). 

The implementation of coupled DE models in GADGET-2 has 
already been described in full detail in BAIO, and since no major 
modifications with respect to that setup are required to carry out the 
time dependent coupling simulations presented in this work, we re- 
fer the interested reader to the more detailed description presented 



Box Size L 

Number of baryonic particles N^ 

Number of CDM particles Nc 

Baryon Mass Af;, 

CDM mass Mc{2 = 0) 

Gravitational Softening ts 



80^-1 Mpc 

512^ 

512^ 

4.82 X \lfh-^ Mq 

2.41 X lO^/i-i M, 

3.5 h-^ kpc 







Table 3. The pai'ameters of the high-resolution hydrodynamical N-body 
simulations discussed in Sec. [5] and|4] 



in BAIO and we only summarize very briefly here the main features 
of our modified algorithm: 

• A table of values of the Hubble function H{a) that is com- 
puted for each model by integrating the system ( I14t is read in by 
the code and the value of H at each timestep is linearly interpolated 
from the table; 

• The mass of CDM particles in the simulation box is cor- 
rected at every timestep according to the mass evolution given by 
Eqn. l llOt . and shown for each model in Fig.[T] 

• The small scale particle-particle (computed with the grav- 
itational oct-tree algorithm implemented in GADGET-2) and the 
large scale particle-mesh gravitational accelerations for CDM par- 
ticles are computed taking into account the additional contribution 
arising from the scalar-field mediated fifth-force, as described by 
Eqn. l l52t : 

• The acceleration of CDM particles computed according to the 
prescription given in the previous point receives at each timestep 
an additional contribution from the extra friction term that appears 
in the expression of H in Eqn. l |52| l. 



As discussed in Sec. 12.31 we want to investigate the nonlin- 
ear evolution of structure formation in the context of a few selected 
cosmological models among all the cosmologies described in Ta- 
ble|2] To this end, we use the modified version of GADGET-2 briefly 
summarized above to run high resolution hydrodynamical N- 
body simulations for the models ACDM, EXP010a2, EXP015a3, 
EXP010e2, EXPOlOeS and EXP015e3. The parameters of the sim- 
ulations are summarized in Table [3] Hydrodynamical forces com- 
puted with the SPH [Smoothed Particle Hydrodynamics) algorithm 
implemented in GADGET-2 are acting on baryonic particles in all 
the simulations, while non-adiabatic processes like e.g. radiative 
cooling, star formation, and feedback mechanisms from Super- 
novae or Active Galactic Nuclei are not included in any of our runs. 



3.2 Initial conditions 

The initial conditions for all the simulations are generated by set- 
ting up a random-phase rea lizatio n of the Eisenstein & Hu power 
spectrum (Eisenstein & Hu 19971) according to the Zel'dovich ap- 
proximation ( Zel 'dovicf^ 1970n . where the normalization amplitude 
of the power spectrum is adjusted to the desired value of erg . In do- 
ing so, we are implicitly assuming that the coupling does not affect 
the shape of the initial matter power spectrum. We are therefore dis- 
carding from our treatment any possible early effect of the coupling 
on the statistical properties of the density field. This is in general 
a reasonable assumption since in all our models the coupling is 
quite small at high redshifts. However, in particular for the expo- 
nential coupling models, such small coupling in the early Universe 
could stiU produce some slig ht tilt in the matter transfer functions 
dMainini & Bonomettdl2007 '). Nevertheless, at the resolution level 
achieved in the present work such tilt was already found to have 
very little impact on the subsequent evolution of structure growth 
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Figure 6. The evolution of the scalar mass rh? 



'/s, 



as a function of the e-folding time N for all the variable coupling models under investigation. 



The three panels refer to the total scalar mass of the models with Pc{<t>) oc a*^! , l3c{<t>) oc Q^, and Pc((j>) ex. e^'^'^'^ , respectively. It is interesting to notice 
that the total effective mass in^ acquires negative values during matter domination in all the three different types of models. However, if the values of order 
~ — 10 found in the former two cases might lead to significant instabilities of the models at large scales, the exponential coupling models show very small 
absolute values of the total scalar mass at all redshifts, withm^ never exceeding 0(1) for all the models considered in our high-resolution N-body simulations. 



as compared to the other effects related to the modified perturba- 
tions evolution (BAIO). Higher resolution studies, however, might 
be significantly affected even by such weak distortions of the small- 
scale power spectrum, and should therefore include a full treatment 
of the early effects of the DE interactions on the matter transfer 
function. 

Initial conditions are therefore generated by perturbing the 
positions of particles from a cartesian grid until the desired den- 
sity field is realized, and then by rescaling the particle displace- 
ments from the grid points with the linear growth factor Dj^ of 
each cosmological model between 2 = and z = 60, which is 
the starting redshift of all our runs. Although this convention on 
the normalization of the power spectrum amplitude is probably the 
most commonly used, other choices are equall y valid and have been 
recently adopted in some related studies (Bal di & PettorinduOld : 
iDe Boni et al]|2010!) . 

All the simulations presented in this work have the same ran- 
dom phases for the realization of the power spectrum in the ini- 
tial conditions, and structures are therefore expected to form in the 
same locations in all the runs. Finally, the velocities of particles are 
set according to linear perturbation theory, by a simple relation with 
the computed initial overdensities, which in Fourier space reads: 



vik. 



if{a)aH5{k, a) 



where the growth rate /(a) is defined as 

dlnD+ 



./(«) 



(60) 



(61) 



dlno 

For ACDM cosmologies, the total growth rate is always well ap- 
proximated by a po wer of th e total matter density {Qc + ^b)~' 
where 7 = 0.55 jPeeblej [l980). This approximation is no 
longer valid for coup l ed DE models in general, as discussed in 
Idj Porto & Amendolal J2008l) and BAIO where alternative phe- 
nomenological fitting formulae for constant coupling models have 
been proposed. Also in the case of time dependent couplings the 
growth rate does not follow the ACDM behavior, and observa- 
tions of the growth of str uctures at different redshifts (see e.g. 
ISean & Tangmatithamll2010. . for a recent analysis of growth data) 
might have the power to detect deviations from the standard value 
of the growth index 7 and put constraints on variable coupling mod- 
els. Such a detection would be a clear indication of some cosmo- 
logical modification of standard gravity. The growth factor D+ and 



the ratio of the growth rate /(a) over the ACDM fitting formula 
Slj{/^, as computed from the numerical integration of Eqs. ( 150151b , 
are shown in the left and in the middle panels of Fig. [T] respec- 
tively, for the models investigated with our high-resolution N-body 
simulations. 

Additionally, in the right panel of Fig.|7]we plot for the same 
set of models the derived evolution of the growth index 7 computed 



7 = 



In/ 



ln(n6 + flc) 



(62) 



It is very interesting to notice the strikingly different behavior 
of 7 for the exponential coupling models EXP010e2, EXP010e3, 
EXP015e3, with respect to ACDM and the phenomenological mod- 
els EXP010a2 and EXP015a3. For the latter, a quasi-constant value 
of gamma very close to the ACDM value of 0.55 is found during 
matter domination, and the curves are practically indistinguishable 
from the fiducial ACDM model. On the contrary, the former models 
show a completely different behavior of 7, that even becomes nega- 
tive for a large period of time and reaches its minimum value in the 
redshift range z ~ 10 — 20. Interestingly, a very similar evolution 
of the growth index parameter 7 to the one shown in the right panel 
of Fig.|7]has been recently reported for some F(R) models of mod- 
ified gravitv JMotohashi et alj201CI : lNarikawa & YamamotduOlOl : 
lApplebv & Welleii2O10n . Given the existence of a well known con- 
formal correspondence between coupled DE models and modified 
gravity theories it would be particularly interesting to investigate 
whether such similar behaviors of the gamma index represent an- 
other element of connection between the two scenarios. 



4 RESULTS OF THE N-BODY SIMULATIONS 

We present and discuss in this section the results of our analy- 
sis of the high-resolution simulations described above. In order to 
study the properties of collapsed structures we first need to iden- 
tify groups and gravitationally bound substructures in our simula- 
tions. To this end , we apply the F riends-of-Friends (FoF) and SUB- 
FIND algorithms JSpringel et afluOOl.) to the particles distributions 
in our simulations outputs. For the FoF computation we use a link- 
ing length of A = 0.2 x d where d is the mean particle spacing. 
As we stressed in Sec. 13.21 all the simulations are started from the 
same random realization of the Eisenstein & Hu power spectrum. 
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Figure 7. Left Panel: The Growth Factor Dj^/a as a function of redshift for the six models selected for N-body simulations. The different colors correspond 
to different models, according to the legend. The red curve, corresponding to the model EXP010a2, is almost completely hidden by the green curve, corre- 
sponding to the model EXP015a3. This shows how these two models are practically indistinguishable from each other. Middle Panel: The evolution of the 
ratio of the growth rate function f(a) = c(lnD_|-/dlna to the ACDM fitting formula f2'^/^ for the same set of models. The black curve, representing 
ACDM, is consistent with a ratio of 1.0, as expected. The scale factor dependent models EXP010a2 and EXP015a3 (red and green lines, respectively), are 
indistinguishable from ACDM for most of the expansion history, but show a progressively faster growth at ^ < 1 — 2. The exponential coupling models, 
instead, have already a sensibly faster growth of perturbations at z as high as 100. Right Panel: The derived 7 index for the growth of density perturbations 
for the same set of models as in the previous two panels. While the scale factor dependent models show a value of 7 consistent with the ACDM value of 0.55 
during matter domination, with only small deviations at low redshift (the slight increase at high z is due to the presence of radiation, and is an expected effect), 
the exponential coupling models show a completely different behavior: the 7 index decreases significantly during most of matter domination, even becoming 
negative at some stage, and reaches a minimum in the redshift range 2 ^ 10 — 20 before growing again towards values comparabl e, but still significantly 
lower, to the ACDM c ase. Similar behaviors have been recently reported for some realizations of F{R) modified gravity theories JMotohashi et alj|201(ll ; 
iNarikawa & Yamamotg 20 10:. Ap pleby & Wellet, 2Qia). 



and structures are therefore expected to form in the same positions 
in all the runs. This will allow us to identify the same objects in all 
the simulations and to directly compare their properties. In BAIO, 
due to the slight tilt in the transfer functions used to set up the initial 
conditions for the simulations of the different coupled DE models, 
it was necessary to apply a selection criterion in order to ensure 
that two structures identified in two different simulations could be 
safely considered as being the same object. This selection crite- 
rion consisted in identifying halos found in different simulations 
as the same object only if the most bound particle of each of them 
lied within the virial radius of the corresponding structure in the 
ACDM simulation. Although in the present work there is no tilt of 
the initial power spectrum shape, we decide nevertheless to apply 
the same selection criterion to our group catalogs, since the differ- 
ent timestepping induced by the different evolution of the effective 
gravitational forces in each run, especially at low redshifts, might 
induce some offsets in the final positions of collapsed structures. 
We restrict this selection procedure to the 300 most massive halos 
in our catalogs, which cover a range of CDM virial masses between 
M200 = 3.46 X 10^2 f^-i j^^ ^„jj j^.^^y^ ^ 3 gg X ioi4 /j-i Mq 

in the ACDM cosmology. 



4.1 Matter power spectrum 

As a first step in the analysis of our N-body runs we compute the 
power spectrum of matter density fluctuations P{k, z) at several 
different redshifts, as well as the separate density power spectra for 
the baryonic and the CDM components, Pj, (k, z) and Pc{k,z). The 
evolution of P{k,z) as a function of scale at different redshifts for 
all our high-resolution simulations is shown in Fig. [8] 

We should remind here that all our simulations are normal- 
ized assuming the same amplitude of the large-scale linear power 
spectrum at the present time. This is done by normalizing with the 
same value of cr8(0) = 0.807 the present amplitude of the Eisen- 
stein & Hu linear power spectrum, and then scaling the values of 



density fluctuations to the starting redshift of the simulations with 
the appropriate growth factor for each cosmological model. As al- 
ready stressed above, we have discarded early effects of the cou- 
pling on the initial power spectrum shape since these have already 
been shown to have a minor impact on the final results at the present 
level of numerical resolution (BAIO) even for larger values of the 
early coupling than those considered here. For higher resolution 
works, however, such effects might become important and should 
be taken into account in the initial conditions of the N-body runs. 

As it can be clearly seen in the first and second panels of 
Fig. |7] the models where the coupling strength evolves as a power 
of the scale factor (EXP010a2, EXP015a3), due to the very fast 
decrease of the coupling towards higher redshifts, present a total 
growth factor which is almost indistinguishable from the ACDM 
one for most of the expansion history of the Universe, and the off- 
set with respect to ACDM is due to the different growth at low 
redshifts (z < 1). On the other hand, the more physically moti- 
vated models with a coupling strength related to the evolution of 
the scalar field (EXP010e2, EXPOlOeS, EXP015e3) show a sen- 
sibly different evolution of the growth factor with respect to the 
ACDM case already at redshifts of the order of z ~ 20. It therefore 
comes as no surprise, as it can be seen in the first three panels of 
Fig.m that at high redshifts the power spectra of the EXP010a2 and 
of the EXP015a3 models are practically indistinguishable from the 
ACDM case, while the exponential coupling models EXP010e2, 
EXP010e3, and EXP015e3, show a slightly lower amplitude of the 
power spectrum at all scales, with a weak enhancement of the effect 
for progressively smaller scales. 

The situation is inverted at lower redshifts, as it can be seen in 
the last three panels of Fig. [8] In fact, due to the constantly faster 
growth of CDM density fluctuations, the exponential coupling 
models show a faster evolution of the power spectrum with respect 
to ACDM, such that between z — 1.0 and z — the gap in the 
power spectrum amplitude with respect to ACDM at large scales is 
progressively reduced and these models catch up the ACDM ampli- 
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Figure 8. The matter power spectrum at different redshifts for the six models studied in our set of N-body simulations. In the bottom part of each panel we 
plot the residuals with respect to ACDM for a more clear understanding of the different evolutions of the power spectrum amplitude at different length scales. 
The different behavior of the two types of coupling evolution, with the exponential coupling models already showing a different amplitude at high z, while the 
scale factor dependent models show a strong increase of power at small scales only at z < 0.5 (as illustrated in detail in the text) is clearly visible in the plots. 
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tude at the present time, while at smaller scales (fc ^ 1.0 /i Mpc ^) 
a residual lack of power is still present at 2 = 0. On the other 
hand, the phenomenological models EXP010a2 and EXP015a3, 
due to the sudden increase of the mutual attraction of CDM par- 
ticles, show a tremendous growth of density fluctuations at scales 
smaller than k ~ 1.0 h Mpc~^ between z = 0.5 and 2 = 0, such 
that at these scales the power spectrum amplitude at low redshifts 
is substantially increased as compared to the ACDM case, with the 
effect becoming progressively larger for smaller scales. The most 
extreme case, given by the model EXP015a3, shows an amplitude 
of the power spectrum at fc ~ 10 h, Mpc~^ that is roughly twice as 
large as for the ACDM cosmology. 

This analysis therefore presents us two quite different behav- 
iors of the two classes of models investigated in our numerical runs. 
On one side, the exponential coupling models have a faster growth 
of the power spectrum amplitude over a rather long period of time, 
and starting with a substantially lower power spectrum amplitude 
at all scales at high redshifts, end up with an amplitude of the large 
scale power spectrum comparable to ACDM at z = 0, nevertheless 
showing a residual lack of power at small scales. On the contrary, 
the phenomenological parametrizations of the coupling evolution 
embedded in models EXP010a2 and EXP015a3 present an almost 
identical growth of linear density perturbations to ACDM during 
most of the expansion history of the Universe, followed by a sud- 
den enhancement of the growth for scales below k ~ 1.0 h Mpc~^ 
between z = 0.5 and 2 = as the coupling strength rapidly in- 
creases towards its present value. This late faster growth turns into 
a strong increase of the power spectrum amplitude at small scales. 
The situation at 2 = 0, shown in the last panel of Fig. [8] gives a 
quite interesting picture, with all the models being practically indis- 
tinguishable from ACDM at the largest scales, as a consequence of 
the common normalization of the linear power spectra at 2 = 0, but 
with a broad range of power spectrum amplitudes at small scales. 

This strikingly different behavior of the evolution of linear 
density perturbations in these two classes of coupled DE mod- 
els could provide ways to observationally distinguish among the 
two scenarios with present and future datasets, and to constrain the 
functional form of the coupling evolution. 



4.2 Baryon-CDM linear and mildly nonlinear bias 

In the context of interacting DE models with a constant coupling 
strength it is now a well established result l lAmendolall200(]| . l2004r) 
that the long-range fifth-force acting between CDM particles in- 
duced by the interaction of the DE scalar field with the CDM fluid 
determines a different growth rate of linear density perturbations 
of CDM with respect to the uncoupled baryonic component. This 
can be clearly understood just by having a look at Eqs. ( 1501 ) and 
( 15 It , where the same spatial distribution of density fluctuations in 
the dominant CDM component sources the CDM perturbation evo- 
lution with a strength Vc times larger than for the baryons. In case 
of a constant coupling /3c (i.e. a constant factor Vc), this different 
growth rate is integrated over the whole expansion history of the 
Universe and induces a sizeable linear bias at all scales between 
the amplitude of density perturbations in the two components. This 
characteristic feature makes the linear bias arising from the cou- 
pling clearly distinguishable, at least in principle, from the hydro- 
dynamical bias arising only at small scales as the Universe becomes 
progressively more structured. The study of such effect within con- 
stant coupling models has been extended to the nonlinear regime 
by numerically following the collapse of a spherical overdensity 



(e.g. by iMainini & Bonometta 2OO67 or by means of N-body sim- 
ulations (Maccioetal. (2004), BAIO), finding that nonlinearities 
enhance the bias between the two components. We want to extend 
here the analysis to the case of time dependent couplings for the 
linear and mildly nonlinear regimes, while the strongly nonlinear 
regime will be studied in Sec. 14.61 

The coupling-induced bias between CDM and baryon density 
fluctuations can be studied by comparing the ratio of the density 
power spectrum amplitudes of the two components in the different 
simulations, defined as: 



R{k,z) 



Pb{k,z) 



(63) 



Pc(k,z) 

In Fig. [9] we plot the evolution of the ratio of R to the ACDM 
case -Racdm for all our simulations as a function of scale, at sev- 
eral different redshifts. Also in this case, there is a clear difference 
between the exponential coupling models and the models where 
the coupling depends on the scale factor. For the former, as it can 
be seen in the upper three panels of Fig.|9l the presence of a non- 
negligible coupling already at high redshifts determines a moderate 
linear bias with a weak dependence on scale already at 2 = 3 and 
down to 2; — 1, with a clear hierarchy corresponding to the val- 
ues of the couplings at high redshifts in the different models, while 
the latter class is practically indistinguishable from ACDM in this 
redshift range. 

At later times (lower three panels of Fig. ^ the situation 
changes very quickly, and the evolution of R{k, 2) becomes much 
more entangled than it has been shown to be for the case of constant 
coupling cosmologies. As the coupling in all the models grows to- 
wards its present value, the scale dependence of R{k, z) starts to 
appear, showing how the bias progressively grows when moving 
from the linear to the mildly nonlinear regime of density fluctua- 
tions. Also in this case, as we showed in the previous section for 
the total power spectrum, the evolution turns out to be much faster 
and with a much stronger scale dependence for the EXP010a2 and 
EXP015a3 models as compared to the exponential coupling mod- 
els. Nevertheless, the exponential coupling model EXP015e3 also 
shows a very quick reduction of the bias ratio between 2 = 1 and 
2 = 0, due to the strong increase of its coupling. 

It is very interesting to notice how all the models that share 
the same final value of the coupling /3o seem to converge, at 2 = 0, 
to very similar values of the bias at small scales, irregardless of the 
type of coupling evolution. The distinction between the two classes 
of models is nevertheless still present at the largest scales of the 
simulations, where the hierarchy of the exponential coupling mod- 
els significantly changes with time, while the scale factor depen- 
dent models roughly retain the initial value of the large scale bias 
ratio i? ~ 1. 



4.3 Halo mass function 

For all of our high-resolution N-body simulations we have com- 
puted the halo mass function based on the groups identified by our 
FoF algorithm. The cumulative mass function for all our runs is 
plotted in Fig. (TO) where each panel refers to a different redshift. 
At 2: = all the mass functions have a similar shape and amplitude 
over the whole mass range covered by our catalog, with a discrep- 
ancy from model to model of the order of ~ 10%, which slightly 
increases at the high-mass end. In particular, it is worth noticing 
here how all the models except the EXP010e2 and the EXP015e3 
show a slightly larger number of halos over the whole mass range 
with respect to ACDM at the present time. 
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Figure 9. The ratio of the bias R{k, z) to the ACDM bias Racdm (k, z) as a function of wavenumber k at different redshifts. As extensively discussed in the 
text, at high redshifts (z > 1.0, upper three panels) the EXP010a2 and EXP015a3 models are practically indistinguishable from ACDM, while the exponential 

coupling models already exhibit an almost constant bias at all scales, with a clear hierarchy corresponding to the different values of the coupling at high z. 
The situation changes dramatically at 2 < 1, where the scale dependence of the bias starts to appear clearly for all the models under investigation, and the 
hierarchy of the exponential coupling models is modified due to the more substantial growth of /3c(0) in the EXP015e3 model. At 2 = 0, all the models with 
the same present value of the coupling /3o interestingly seem to end up with comparable values of the bias at the smallest scales available to our analysis. 



At higher redshifts, instead, the models under investigation 
show very different mass function evolutions. On one side, the 
EXP010a2 and EXP015a3 models show very little differences 
from ACDM. On the other side, the exponential coupling models 
have always a significantly lower number of halos with respect to 
ACDM at intermediate and high masses, while a slight excess of 
low mass halos is clearly visible for these models at z ~ 1 — 2. 
This behavior is due to the fact that in the exponential coupling 
models structure formation is starting later than in ACDM as a con- 
sequence of having normalized all the cosmologies to the same as 
at the present time. This shows how in these cosmologies halos 
of any given mass tend to form later than in ACDM, and how at 
2: ~ 1 — 2 most of the small halos did not have the time yet to 
merge and form larger and more massive structures. However, as it 
can be clearly seen by the evolution of the mass function with red- 
shift, the gap in the number of large halos between the exponential 
coupling models and ACDM is progressively reduced as time goes 
by due to the higher growth rates of the coupled cosmologies, as the 
increase of the CDM mutual attraction speeds up the aggregation 
of small objects into larger structures. 

We have also computed the multiplicity function (defined as 
[M^/p]-dn(< M)/dM) for all of our models, which is plotted 
in Fig . [TTI where we also plot fo r comparison the lSheth & TormenI 
( 1 19991) and Ijenkins et al.l ( 1200 1[) fitting formulae evaluated at dif- 
ferent redshifts using the appropriate growth factor for each model, 
and with the standard value of the extrapolated linear density con- 
trast at collapse Sc — 1.686. We find that the usual mass function 
formalism reproduces fairly well the distribution of our simulated 



halos up tp 2 = 3 in most of the investigated models, extending 
the validity of recent findings for early dark energy cosmolo gies 
JGrossi & Springelll2009l : lFrancis et al.ll2008l : |Pace et al.ll2010h and 
for constant coupling interacting DE models (BAIO), to the case of 
variable couplings. However, some sizeable discrepancy between 
the predicted multiplicity function and the outcomes of our sim- 
ulations appears for the EXP015a3 model at z — 0, where both 
the mass functions fitting formulae systematically underestimate 
the simulated distribution. This behavior might be related to the 
sudden and strong increase of halo masses at z ~ 0.5 — in 
this model, which would explain why the offset appears only at 
z — 0. However, a careful analysis of the spherical collapse for- 
malism in the context of variable coupling models of interacting 
DE, extending t o these scen arios the analysis recently carried out 
by Winterg erst & Pettoring (2010) for constant couplings, will be 
necessary in order to fully understand the reasons of this discrep- 
ancy. We defer such analysis to future work. 



4.4 Halo density profiles 

Cosmological simulations of structure formation have consistently 
shown that the density profiles of dynamically relaxed CDM halos 
have a universal shape that can be accurately fitted for any hal o 
mass by the NEW fitting formula ( iNavarro et al.1 19911 1996l.ll997h : 



pcrit 



(t) 



1 + 



(64) 
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Figure 10. The mass functions of all the six models studied with our set of N-body simulations, plotted at different redshifts. At the present time all the models 
show similar mass functions, with differences of the order of '~ 10% with respect to the ACDM case, slightly increasing at the high-mass end. At higher 
redshift the differences among the models become more pronounced, and the exponential coupling models EXP010e2, EXP010e3, and EXP015e3 show a 
slight excess of low mass objects and a considerable lack of high mass objects at redhifts in the range z ~ 1 — 2, while the remaining models do not show 
very significant differences with respect to ACDM. This picture is consistent, as explained in the text, with the fact that in the exponential coupling models 
structure formation starts later due to the lower initial amplitude of the power spectrum. The gap is then reduced at redshifts between 1 and as a consequence 
of the strong increase of the growth factor in these models. 



where pcrit = 3H^M^ is the critical density of the Universe, 5c is 
the characteristic halo density contrast, and r^ is the halo scale ra- 
dius. However, several astrophysical observations of mass density 
profiles of dwarf Low Surface Brightness (LSB) galaxies JMoord 
1 1994 lFlores&Prim ack"l994': 'Simon et al."'20 03|), of Milky- Way 
type spiral galaxies 1 Navarro & Steinmetz^200Q : ISalucci & Burkerd 



l2000l:rSaluccii200ll. Binnev & Eyansi200lh, oreven o f large galaxy 
clusters iSand et alj2002ll2004lNewman et alj2009l) . have shown 
at different levels that these objects have shallower density pro- 
files than predicted by the theoretical universal NFW shape. This 
tension between simulations and observations is often referred to 
as the "cusp-core" problem. Several attempts have been made in 
order to solve this discrepancy by invoking backreaction mech- 
anisms of the baryonic compo nent on the CDM density pro- 
files (see e.g. iDuffv et al.l 120100 or by different flavors of warrn 
dark matter (WDM) JAvila-Reese et al.ll200ll : IStrigari et alj|2007l : 
Ide Narav et alj2()09l:lde Vega et al.l20ia) and Self Interacting Dark 
Matter (SIMD) JSpergel & Steinhardtll2000l : IWandelt et alJl2000l : 



iDave et al.l200lh . Here we want to investigate another independent 
possibility, namely the fact that a DE-CDM interaction could play 
a role in alleviating the "cusp-core" tension. 

One of the main results found in BAIO for the case of constant 
couplings consisted in the discovery that the nonlinear dynamics of 
interacting DE cosmologies determines a systematic reduction of 
the inner overdensity of CDM halos with respect to ACDM, with 
the effect growing for increasing coupling. This result, in stark con- 
trast with previous works, was the first evidence of how interacting 
DE models could produce shallower density profiles and less con- 
centrated halos, thereby providing a possible solution to the "cusp- 
core" problem. Nevertheless, the present observational constraints 
on constant coupling models (BE08,LV09) put tight bounds on the 
maximum allowed amplitude of this effect, which turns out to be 
not strong enough to fully address the problem. In particular, it is 
worth reminding here that the largest coupling considered in BAIO 
(/3c = 0.25), which was found to determine a reduction of the inner 
overdensity by ^ 20%, is already observationally ruled out even by 
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Figure 11. The multiplicity functions for the six models investigated with N-body simulations. In each panel, the differently-coloured sets of data points, with 
the relative error bars, are the multiplicity functions evaluated in equally spaced logarithmic mass bins at four different redshifts. The dashed and dot-dashed 
lines represent the predictions for the multiplicity functions according to the Sheth & Tormen 1 1999) and Jenkins et al. (2001) fitting formulae, respectively. 
The plots show a fairly good agreement between the simulated multiplicity function and the theoretical predicitions, with the sole exception of the EXP015a3 
model at z = {right upper panel), where both the fitting functions underestimate the simulated distribution. 



the weaker bounds derived by LV09 ( | /3c | < 0.17) for the case of a 
large average neutrino mass. 

It is therefore very interesting to investigate whether variable 
coupling models - where the coupling is relatively small for a large 
fraction of the expansion history of the Universe - could produce 
stronger effects on the density profiles without running into conflict 
with the present observational bounds, as it is the case for all the 
models selected for N-body simulations in the present work, which 
are at least consistent (according to the selection criterion described 
in Sec. |2.3l l with the constraints derived by LV09. 

To this end, we have computed the spherically averaged total 
matter (baryons and CDM) density profile as a function of radius 
around the halo center (defined as the position of the particle with 
the minimum gravitational potential within the group) for all the 
halos in our group catalogs that can be safely identified as the same 
object arising in the different simulations. In Fig.[T2]we plot the to- 
tal matter density profiles for four halos of different mass in the six 
cosmological models under investigation. The left upper panel of 
Fig. [12] shows the highest mass halo in our sample, while the right 
lower panel show the lowest mass one. Interestingly, we find that 
not all the models show a decrease of the inner halo overdensity as 
it was found for the case of constant couplings, where the ACDM 
model always showed a larger overdensity than any coupled DE 
cosmology. On the contrary, some of the models (the ones where 
the coupling depends on a power of the scale factor), are even found 
to have significantly larger values of the central overdensity with 
respect to ACDM, therefore showing an opposite trend to the one 
required in order to address the "cusp-core" problem. On the other 
hand, the exponential coupling models show the expected system- 
atic lowering of the inner density, although the effect is clearly not 



yet strong enough in these models to produce a cored profile. Once 
more, we have found a very different behavior of the two classes of 
variable couplings under investigation, and the interesting situation 
depicted in Fig.[T2]deserves to be extensively discussed. 

The reason for these strikingly different evolutions of the non- 
linear dynamics of CDM halos within the two different classes of 
models can be understood by having a look at Fig.|4] In BAIO it was 
clearly shown how the reduction of the inner overdensity of halos 
in constant coupling models (and the consequent reduction of halo 
concentrations, that will be discussed in the next section) was pri- 
marily due to the effect of the friction term IjScXVc in Eqn.[52] In 
particular, it was shown that the energy gained by the collapsed sys- 
tems due to the extra acceleration of CDM particles induced by the 
friction term could move the systems out of their virial equilibrium 
and determine a slight expansion of the halos. One of the conse- 
quences of such slight adiabatic expansion is the transfer of mass 
from the center of the halos towards the outskirts, thereby deter- 
mining a reduction of the inner overdensity. It is therefore clear, by 
looking at Fig.|4] why this mechanism does not produce the same 
effects for the EXP010a2 and EXP015a3 models as it does for the 
constant coupling models studied in BAIO: although having large 
values of the coupling at low redshift, these models feature a quite 
small friction term due to the absence of a "^MDE". This deter- 
mines a fast decay with redshift of the scalar field kinetic energy 
X and a consequent strong suppression of the extra friction term 
during most of the cosmological evolution, as discussed in detail 
in Sec. 12.41 On the contrary, in the exponential coupling models 
EXP010e2, EXPOlOeS, and EXP015e3, the presence of what we 
called a "Growing (jiMDE" is able to sustain a slower decay of 
the scalar field kinetic energy with redshift, and the friction term 
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Figure 12. Total matter (CDM + baryons) density profiles for four halos of different masses in the six models investigated with our set of N-body simulations. 
The left upper panel represents the most massive halo in our sample, while the right lower panel shows the lowest mass one. The vertical dot-dashed line 
indicates the location of r2oo for the ACDM halo. In all the four plots it is clear how the variable coupling models investigated in the present work, contrarily 
to what happens for constant couplings, do not always determine a decrease of the inner overdensity of halos. In particular, as discussed in full detail in the 
text, the models where a "c/iMDE" or a "Growing t/iMDE" is absent show a significant increase of the inner overdensity. 



is therefore less suppressed and still capable of inducing the ex- 
pansion of CDM halos and the consequent reduction of the inner 
overdensity. 

We have therefore shown here that the presence of a "0MDE" 
or of a "Growing (/>MDE" is not only a desirable feature of an in- 
teracting DE model due to its capability of easing the "fine tuning 
problem", but is also an essential ingredient in determining the type 
of impact that the interaction can have on the nonlinear dynamics 
of coupled matter particles at small scales. 

Even though we have now explained why the absence of a 
"Growing ^MDE" suppresses the efficiency of the coupling in low- 
ering the density profiles of collapsed structures, we still need to 
explain why some of our models even show a significant increase 
of the inner overdensity of halos over the whole mass range of our 
catalog with respect to the uncoupled ACDM case. The explana- 
tion of this new effect seems also clear if one considers the virial 
equilibrium of a collapsed halo. Along with the gain of energy due 
to the extra friction term, there are two other effects that can mod- 
ify the virial state of a system within interacting DE cosmologies. 



On one side, the variation of the mass of CDM particles described 
by Eqn. [TO] determines a progressive increase of the gravitational 
potential energy of the system. The effect of this mass loss would 
also be a slight expansion of the halos. However, as it was shown 
in BAIO, this effect is expected to be rather small for models with 
an overall variation of the CDM particle mass comparable or even 
larger than for the specific models under investigation here. On 
the other side, the fast growth of the effective gravitational con- 
stant acting between CDM particles G{(l)) = Tc{(t))G determines a 
corresponding fast decrease of the gravitational potential energy of 
the system that balances and exceeds the small increase due to the 
mass variation. This latter mechanism is not present at all in con- 
stant coupling models, but plays a crucial role here where the factor 
Vc{(j)) can grow by ~30 to 75 % during cosmic evolution. The total 
combined effect is therefore a net decrease of the total energy, and 
a consequent contraction of the halos forming in the models where 
the friction term is suppressed and cannot provide to the systems 
the energy increase necessary to counteract this mechanism. Con- 
sistently with this picture, halos are found to be more overdense in 
the EXP010a2 and EXP015a3 models as compared to ACDM, and 
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Figure 13. NFW fit (grey, dot-dashed fines) to the density profiles (sofid fines) of the most massive {upper three panels) and the least massive {lower three 
panels) halos in our sample in the ACDM {left panels), EXP010e2 {middle panels) and EXP015a3 {right panels) cosmological models. Each plot reports the 
halo mass A/2ao, the halo radius r2oo, the scale radius r^, the concentration c, and the x^ of the NEW fit. A clear trend of increase of the scale radius in the 
exponential coupling models, and of decrease in the scale factor models appears at both mass scales. 
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Group 278 Group 278 
r^(h-ikpc) ^^(ACDM) 



184.811 
151.660 
99.528 
295.849 
213.840 



1.0 
0.821 
0.539 
1.601 
1.157 
1.211 



112.022 
85.959 
61.153 
174.661 
127.927 
130.655 



1.0 
0.767 
0.546 
1.559 
1.142 
1.166 



51.244 
44.900 
32.913 
75.177 
59.603 
62.153 



1.0 
0.876 
0.642 
1.467 
1.163 
1.213 



Table 4. Evolution of the scale radius Vs for the four halos shown in Fig. [12] with respect to the corresponding ACDM value. Contrarily to what happens 
for constant coupling models the scale radius can either increase or decrease with respect to ACDM in variable coupling scenarios according to the type of 
coupling evolution: models that do not present a "Growing r/iMDE" phase feature a contraction of collapsed halos and a consequent decrease of their scale 
radius up to ~ 45%, while the exponential coupling models show an increase of the halo scale radius up to ~ 60%, a much more significant effect than for 
constant coupfing models. 



less overdense in the other exponential coupling models, where the 
friction term balances - and in some cases exceeds - the decrease 
of gravitational potential energy. 

It is therefore clear that a significant growth in time of the ef- 
fective gravitational constant in any cosmological modification of 
newtonian dynamics will result in an increase of the central den- 
sity of CDM halos unless there are other mechanisms - as e.g. a 
large extra friction term in the case of the interaction of CDM 
with the DE - able to balance the decrease of gravitational po- 
tential energy. Hence, it is natural to expect that also other cos- 
mological models that introduce an effective modification of the 
gravitational interaction due to a long or a sho rt-range fifth-force 
(as e.g. the recently proposed ReBEL scenario: iGubser & Peeblea 



2004lFarrar & PeebleslZOoi i lParrar & Rosenl2007l : lHelIwing et all 



201(J) will produce larger values of the halo inner overdensity and 
more cuspy profiles with respect to ACDM if the strength of the 
additional fifth-force grows in time during the epoch of structure 
formation. 



The last point that we still need to address in this section is 
whether or not the time variation of the coupling allows to have a 
more significant reduction of the inner overdensities of collapsed 
objects with respect to the constant coupling case, without running 
into conflict with observational constraints. We have already shown 
that none of the models investigated in this work actually presents 
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Figure 14. The ratio of the spherically averaged matter density around the 
center of a halo to the density of the same halo in the ACDM cosmology, 
averaged over all the halos in our sample and plotted as a function of frac- 
tional radius with respect to r20o ■ The plot clearly shows how to a reduced 
inner density corresponds a higher density in the outer parts of the halos, 
due to the transfer of mass from the central regions to the outskirts. The 
plot also shows how in some cases (e.g. for EXP010e2 model among the 
cosmologies considered here) a time dependent coupling can produce com- 
parable or even stronger reductions of the inner overdensity of halos with 
respect to constant coupling models, without determining the same level of 
tension with present observational constraints on the background expansion 
of the Universe. 



cored density profiles. This leads us to the conclusion that a com- 
plete solution of the "cusp-core" problem cannot be achieved by the 
simple classes of coupling evolution considered here. Nevertheless, 
at least the EXP010e2 model is found to determine a significant 
lowering of the halo density profiles in the inner regions. In order 
to quantify this effect, and to compare it with the constant coupling 
case, we have computed for our sample the average radial density 
ratio with respect to ACDM, defined as: 

PM{r) 



, w , (65) 

where the average is taken over all the halos in the sample. This 
quantity is plotted in Fig. [14] for all the models under investiga- 
tion and for comparison also for the RP5 model studied in BAIO 
as a function of the fractional radius r/r2oo, where r2oo is the ra- 
dius enclosing a mean overdensity 200 times larger than the critical 
density pcrit- By looking at Fig. [14] it is first of all interesting to 
notice how all the models that show an average lower density in 
the inner regions of CDM halos with respect to ACDM, show cor- 
respondingly larger values of the density in the outskirts. This be- 
havior witnesses the fact that the lower inner overdensity is due to 
the transfer of mass form the core of the halos to the outer regions, 
as it was demonstrated in BAIO for the case of constant coupling. 
The opposite applies to the scale factor dependent models where 
the contraction of the halos brings mass from the outer regions to 
the center. 

Finally, we notice how indeed the time dependence of the 
coupling can produce a comparable and in some cases a stronger 
reduction of the halo inner density with respect to constant 
coupling models, without determining the same impact on the 
overall background evolution of the Universe. In particular, the 
EXP0I0e2 model determines an average decrease of the inner 



overdensity with respect to ACDM roughly 1.5 times larger 
than the RP5 constant coupling model. For the other exponential 
coupling models EXPOlOeS and EXP015e3 the balance between 
the friction term and the decrease of gravitational potential energy 
is less favorable, and they are found to produce weaker effects with 
respect to RP5, although still having a much lower impact on the 
background expansion, which makes them still preferable over the 
constant coupling scenario. 

Despite all these significant effects on the inner overdensity 
of CDM halos, we find that the density profiles can be still fitted 
extremely well with an NEW shape in all the models. This striking 
result confirms also for the case of time dependent couplings, as 
it was found for constant couplings, that the shape of the density 
profiles is not significantly modified by the additional physical pro- 
cesses related to the DE-CDM interaction, and that the change of 
inner overdensity is related to a change in the location of the scale 
radius r^ in the different cosmologies. This can be clearly seen in 
Fig. [13] where we plot for the largest and the smallest halos in our 
sample (upper and lower three panels, respectively) the density pro- 
file and the best-fit NEW function for ACDM and for the two most 
extreme models in each of the two different classes of coupling evo- 
lution, the EXP010e2 and EXP0I5a3 models. As Fig.[T3]shows, the 
NEW fit to the simulated halos is equally good (the x^ for the fit 
is indicated in each figure) at the two mass extremes of our sam- 
ple in all the models. Consistently with our interpretation, there is 
a clear trend of the scale radius rs which always increases in the 
EXP010e2 and decreases in the EXP0I5a3 models. This effect is 
also clear from Table|4] where we list for the four halos considered 
in Eig.[T2]the scale radius in all our simulated cosmologies, and its 
ratio to the ACDM value. Although some exceptions to this general 
trend can be found in our halo sample, Table|4]shows the tendency 
of the vast majority of the halos, with a decrease of Ts up to ~ 45% 
in the scale factor dependent models, and an increase up to ~ 60% 
in the exponential coupling models, with respect to ACDM. 



4.5 Halo concentrations 

As anticipated in the previous section, we also compute for all the 
objects in our sample the halo concentrations with two independent 
methods. First, we compute the concentration of a halo as: 

'"200 .£^^ 

c = (66) 

based on a NEW fit of the density profiles . Then, we use the inde- 
pendent method devised bv iSpringel et alj OOOSh to compute halo 
concentrations according to the equation: 



200 



3 ln(l + c) - c/(l + c) 
with 5v defined as: 

ymax 



7.213 5v 



Sv ^2 



HoT-r, 



(67) 



(68) 



where Vmax and Vmax are the maximum rotational velocity of the 
halo and the radius at which this velocity peak is located, respec- 
tively. Following the notation used in BAIO we will denote the con- 
centrations computed with this method as c* . 

The results of these two independent methods to compute con- 
centrations are shown in Fig. [15] where concentrations are plotted 
as a function of the halo virial mass A/200 for all the high-resolution 
simulations carried out in the present work. Despite the different 
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Figure 15. Evolution of the mean halo concentrations as a function of mass for the 300 most massive halos in our group catalogs, in the different cosmological 
models under investigation. Th e concentrations have been computed by a direct fit of the density profiles with an NFW shape (left panel) or by using the 
independent method devised by Springel et al. (2008) and described in the text (right panel). The filled circles represent the mean halo concentrations in each 
of the four mass bins in which the mass range of our group catalogs has been divided. The con'espondingly coloured dotted lines indicate the spread of 68% 
of the halos in each mass bin. Both plots show the same trend of concentrations in the different cosmological models, with the exponential coupling models 
being less concentrated than ACDM while the remaining models show significantly larger values of concentrations over the whole mass range. 



methods used, the two panels of Fig. [15] show the same relative 
trend for the concentrations in the different cosmological models. 
Consistently with what found for the halo density profiles in the 
previous section, the halos formed within one of the scale factor 
dependent coupling models are found to have significantly higher 
concentrations than ACDM halos at all masses. Once again, we see 
here that in case of time dependent couplings the ACDM model 
no longer represents an extreme for the range of models investi- 
gated, as it happens for constant couplings. On the contrary, devi- 
ations from the standard cosmological model are possible in both 
directions for different types of coupling evolution. The exponen- 
tial coupling models, as expected, show lower halo concentrations 
with respect to ACDM, and the hierarchy of models follows the 
same order already shown for the iimer overdensity of CDM halos 
in the previous section. 



4.6 Nonlinear bias and halo baryon fraction 

Another source of tension between the predictions of the ACDM 
paradigm and a number of astrophysical observations concerns the 
baryonic budget of large galaxy clusters. The baryon fraction con- 
tained in massive clusters of galaxies is expected to be a fair sample 
of the background cosmological baryon fraction. However, several 
observ ational estimations o f the baryonic content of X-ray clusters 
(as e.g . lAUen et alj (l2004l); lyikhlinin et all (l2006h: IL aRogue e£ai] 
( 120061) ; I Afshordi et al.l ( I2OO7I) . but see also lGiodini eTalJ ( I2009D for 
a recent opposite claim) seem to indicate that these objects have a 
lower content of baryons as compared to the background baryon 
fraction estimated fr om cosmological observations as e.g. CMB 
JKomatsu eraLll2010h . 

As we have already shown in Sec. 14.21 one of the characte- 
ristic features of interacting DE cosmologies is the gravitational 
bias that develops between the amplitude of density fluctuations of 
baryons and CDM due to the different effective gravitational dy- 
namic equations (see Eqs. |52|5Jt that coupled and uncoupled par- 
ticles obey. As a consequence of this different evolution, the bary- 
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Figure 16. Nonlinear bias between the overdensity of baryonic and CDM 
particles as a function of distance from the center of one of the halos shown 
in Fig. 1 121 in all the six different cosmologies considered in our set of N- 
body simulations. A similar trend is found for all the other halos in our 
sample. The enhancement of the bias in the innermost regions of the halo 
is evident for all the cosmologies, and becomes progressively stronger for 
larger values of the present coupling strength /3o. The vertical dot-dashed 
line represents the location of the virial radius r2oo of the halo in the ACDM 
cosmology. 



onic fraction of any overdense region of the Universe is no longer 
expected to match the background cosmological value, even in the 
linear regime, as we showed in Sec 14. 21 

We want to extend here the analysis of this effect to the non- 
linear regime of structure formation, and give an estimate of how 
the halo baryon fraction at z = can be affected by time depen- 
dent couplings in the dark sector. To do so, we first compare the 
ratio of baryon to CDM overdensity as a function of radius around 
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different models (represented by the open diamonds) tliat it is not 
particularly rare in these latter models to find halos with relative 
baryon fractions of the order of ~ 0.7, in particular at intermediate 
and small masses. 

We stress again here that our simulations include hydrody- 
namical forces for the baryonic particles, but do not include other 
non-adiabatic processes like e.g. radiative cooling, star formation, 
and feedback. Therefore, we do not expect our predictions for the 
baryonic fraction of halos to be directly comparable with observa- 
tions; nevertheless, it is clear that the strong reduction of the bary- 
onic content in collapsed objects shown in Fig.[T7]would still be in 
place also in the presence of such non-adiabatic processes, which 
would therefore operate on a significantly smaller reservoir of bary- 
onic mass. These results therefore indicate that a time dependent 
interaction between DE and CDM might be considered as one of 
the possible explanations for the observed low baryon fraction of 
galaxy clusters. 



Figure 17. Relative baryon fraction Yf, within the virial radius r200 for the 
300 most massive halos in our group catalogs plotted as a function of virial 
mass Af2oo. The open diamonds represent the relative baryon fraction of 
individual halos for a uniform random sampling of our catalogs, while the 
filled circles show the mean relative baryon fraction in each of the five mass 
bins in which the full available mass range has been subdivided. Dotted 
lines indicate the spread around the mean of 68% of the halos for each 
cosmological model. A strong reduction of the relative baryon fraction is 
clearly visible for all the coupled DE models as compared to ACDM, with 
the total effect being approximately proportional to the different values of 
the coupling at the present time /3o . 



the center of a halo, defined as: 



B{< 



Pbi< r) ~ pb 



Pc 



(69) 



pb pc(< r) — /9c 

for all the halos in our catalog, as it was done in iMaccio et al.l 
1 120041) and BAIO for the constant coupling case. In Fig.[T6]we plot 
the ratio B{< r) as a function of radius for one of the four halos 
already shown in Fig. [12] The strong enhancement of the bias when 
moving towards the center of the halo follows a very similar behav- 
ior to the case of constant coupling models studied in lMaccio et al.l 
( 120041) and BAIO, although the amplitude of the effect is gener- 
ally larger, consistently with the higher values of the coupling at 
low redshift in the models considered here as compared to previous 
works. 

We then compute the relative baryon fraction Yt for all of our 
halos, defined as: 



Yb 



Mb{< r2oo) 






(70) 



Mtot{< r2oo) 

which is plotted in Fig.[T7]as a function of halo virial mass. While 
the ACDM baryon fraction is in full agreement with previous find- 
ings for the ACDM model (Et tori et alJ l2006':'Gottl oeber & Yepej 
|2007|) . the coupled DE models show as expected a significant re- 
duction of the relative baryon fraction over the full mass range 
covered by our catalogs of simulated halos, with average baryon 
fractions that decrease down to a value of Yt ~ 0.76 — 0.81 for 
the most extreme cases represented by the models EXP015e3 and 
EXP015a3. These values are somewhat lower, as expected, than for 
the constant coupling scenarios investigated in previous works. It is 
also interesting to notice, as it can be seen from the spread of the 
68% of the halos (indicated in Fig. [T7]by the dotted curves above 
and below the mean) or from the location of single halos in the 



5 CONCLUSIONS 

In the context of interacting DE models we have studied a few gen- 
eral classes of time evolution of the interaction strength between 
DE and CDM, generalizing the widely studied case of constant cou- 
plings to the more natural scenario of a variable coupling. 

Following the idea - already discussed in previous works - 
that a large value of the coupling could leave distinctive features in 
the properties of observable structures and even possibly alleviate 
the tensions between the ACDM cosmological model and astro- 
physical observations at small scales, we have designed a few gen- 
eral forms of coupling functions /3c (<^) that grow in time, thereby 
having a significantly weaker impact on the overall background ex- 
pansion of the Universe as compared to constant coupling models 
with the same interaction strength at 2 = 0. 

In particular, we have investigated three classes of time evolu- 
tion of the coupling, where the interaction strength is proportional 
either to a power of the cosmological scale factor a{t), or to the 
fractional DE density Sl^, or to an exponential function of the DE 
scalar field (jf>. The first two classes are purely phenomenological 
parametrizations of the time evolution of the coupling, while the 
latter one represents a more physical situation where the interaction 
depends on the dynamical evolution of the scalar field. 

We have performed a complete numerical analysis of the back- 
ground evolution for these three different types of coupling func- 
tions by solving the full system of coupled dynamic equations in the 
presence of a variable coupling, generalizing previous works. Even 
in the absence of analytic solutions for variable coupling models, 
our numerical integrations allow to identify the main background 
features of these cosmologies. 

More specifically, we have shown that the first two phe- 
nomenological parametrizations of the coupling evolution 
mentioned above, due to the very fast decrease of the coupling 
with increasing redshift, do not present the so called "0MDE" 
scaling solution typical of constant coupling models, and for what 
concerns their background evolution are practically indistinguish- 
able from ACDM, thereby suffering of the same level of fine 
tuning of the cosmological constant. On the contrary, the more 
physical form of a coupling that depends on the evolution of the 
scalar field shows a background evolution with an intermediate 
behavior between ACDM and a standard "<^MDE'" phase which 
can be still well reproduced by the usual analytic solution for 
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the fractional DE density Q.^ during the "(/>MDE" phase, once 
generalized to the case of growing couplings. We have therefore 
called this intermediate type of background evolution a "Growing 
(^MDE" phase. 

We have then studied the evolution of linear perturbations 
within variable coupling models, pointing out the main differences 
arising in the perturbations equations due to the time dependence 
of the coupling. In particular, we have shown how in general a 
growing coupling function could induce instabilities in the growth 
of scalar perturbations at large scales, due to the presence of a 
negative effective mass term in the linear scalar field equation, and 
we have discussed under which conditions these instabilities can be 
avoided. We have then numerically computed the growth factor of 
matter density perturbations at subhorizon scales for a few selected 
models. Among these, the computed growth factors again clearly 
shows two distinct classes: on one side the phenomenological 
couplings present no significant differences in the growth of 
density perturbations with respect to ACDM, except at very low 
redshifts, while on the other side the exponential coupling models 
have a faster growth of density fluctuations during most of the 
expansion history of the Universe. 

The main focus of the present paper is on the effects of 
variable couplings on nonlinear structure formation. Exploiting 
the implementation of coupled DE cosmologies into the N-body 
code GADGET-2 developed for previous works, we have run high- 
resolution hydrodynamical N-body simulations for some selected 
cosmological models belonging to different classes of coupling 
evolution. For all these cosmologies initial conditions have been 
generated discarding possible early effects of the coupling based 
on the consideration that these effects were shown by previous 
studies to have a minor impact on the final properties of nonlinear 
structures at the present level of numerical resolution. We have 
also chosen to normalize all the cosmologies to the same amplitude 
of the large-scale power at the present time. Although this is a 
coiranon choice, other conventions are equally valid and could 
lead to different predictions for the same cosmological models. 

We have shown that the power spectrum of matter density 
fluctuations evolves in a strikingly different way in the two differ- 
ent types of models with respect to ACDM. In particular, the phe- 
nomenological couplings have a very similar evolution to ACDM 
until z ~ 0.5, followed by a very fast growth of the power spec- 
trum amplitude at scales below k ~ 1.0 /i Mpc^^. On the contrary, 
the exponential models have a lower amplitude than ACDM at all 
scales during most of the cosmic evolution, and catch up with the 
ACDM power spectrum only at large scales at 2 = 0, while at small 
scales they still show a significant lack of power also at the present 
time. 

We have then confirmed that also for variable couplings the 
scalar fifth-force acting only between CDM particles induces a bias 
in the amplitude of density fluctuations in baryons and CDM at all 
scales, as it happens for constant coupling models. The bias shows a 
clear scale dependence that develops very quickly as the couplings 
approaches their large values at low redshift, and the enhancement 
of this effect when moving from the linear regime of very large 
scales to smaller and progressively more nonlinear scales is found 
to be stronger than for the constant coupling models studied in pre- 
vious works. 

This bias can be detected also in the very nonlinear regime 
characterizing the inner parts of collapsed objects, and has an im- 



pact on the total amount of baryons contained in massive halos that 
could therefore influence the determination of the baryon fraction 
from cluster measurements. We have therefore computed the evo- 
lution of the average baryon fraction within the virial radius r2oo 
for all the halos arising in the simulations of the different cosmo- 
logical models, finding a generally stronger reduction of the baryon 
fraction as compared to constant couplings, with a decrease up to 
~ 14 - 16% with respect to ACDM. 

We have computed the mass functions at different redshifts 
for all of our selected models, showing how at 2 = all the cos- 
mologies have similar shapes and amplitudes of the mass function, 
with relative differences of the order of ~ 10%. However, at higher 
redshifts the mass functions of the exponential coupling models 
show a clear excess of small halos and a strong lack of large halos 
with respect to ACDM, consistently with the later onset of structure 
formation in these models. We have also computed the multiplicity 
function for all the mode ls and compared i t with t he t heoretical pre- 
diction s according to the lSheth & TormenI ( 119990 and ljenkins et al.l 
( 1200 II) fitting formulae, evaluated with the appropriate growth fac- 
tor for all the models, and with the standard value of the extrapo- 
lated linear overdenisty at collapse. We found good agreement be- 
tween the multiplicity functions in our simulated cosmologies and 
the theoretical fitting functions, with the only exception of one of 
the scale factor dependent models at 2 = 0. This discrepancy might 
suggest the need to reconsider the spherical collapse formalism in 
the presence of strongly variable couplings, will be investigated in 
future works. Nevertheless, the usual mass function fitting formu- 
lae were found to be fairly accurate also for most of our variable 
coupling models, generalizing previous results. 

Finally, we have investigated the effects of variable couplings 
on the halo density profiles. As for the case of constant coupling 
models, we find that they are still remarkably well fit in all the 
different cosmologies by the NEW formula. However - in con- 
trast with what happens within constant coupling models - we also 
find that variable coupling cosmologies do not always show a de- 
crease of the inner overdensity of halos with respect to the stan- 
dard ACDM case, but present opposite trends for the two differ- 
ent classes of coupling functions. The more realistic and physically 
motivated exponential coupling models show a significant decrease 
of the inner overdensity of halos with respect to ACDM, while the 
phenomenological models show on the contrary a clear increase of 
the density in the central regions. This strikingly different behavior 
can be explained by considering which are the main physical mech- 
anisms that can account for a modification of the equilibrium state 
of a collapsed halo in the context of our variable coupling mod- 
els. As it was shown before for the case of constant couplings, the 
mass decrease and the extra friction term in the equation of motion 
of CDM particles can only induce an increase of the total energy 
of a virialized system, which therefore restores its virial equilib- 
rium by slightly expanding. This effect is the source of the lower 
densities in the cores of CDM halos in the presence of constant 
couplings. However, if the coupling is changing in time, there is 
an additional mechanism coming into play: the total potential en- 
ergy of the system decreases due to the increase of the effective 
gravitational constant as a consequence of the growing scalar fifth- 
force. In our models, the effective gravitational constant grows by 
~30 - 75 % during the whole cosmic evolution, which determines 
a corresponding decrease of the gravitational potential energy of 
collapsed systems. This decrease of total energy then determines a 
contraction of the halos. 

Interstingly, we have shown that the two opposite behaviors 
found for the inner overdensity of nonlinear structures are deter- 
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mined by the background evolution of the coupled DE-CDM sys- 
tem: in the models that do not present a "Growing (^MDE" phase 
the friction term, which is the main driver of the expansion of ha- 
los in constant coupling models, is strongly suppressed and cannot 
counteract the contraction induced by the strong increase of the ef- 
fective gravitational constant. On the contrary, the exponential cou- 
pling models, due to the presence of a "Growing ^MDE", have a 
still efficient friction term that balances and sometimes overcomes 
the effect of the potential energy decrease, thereby determining an 
overall expansion of the halos. 

We have therefore shown how the nonlinear behavior of matter 
particles at small scales, in the context of coupled DE models with 
time dependent couplings, is directly influenced by the cosmologi- 
cal background evolution of the scalar field, and how the presence 
of a "Growing ^MDE" phase is essential to determine whether ha- 
los will contract or expand in these cosmologies. 

These considerations apply also to halo concentrations, which 
are found to be higher with respect to ACDM in the phenomeno- 
logical models due to the absence of a "Growing <jf>MDE", and 
lower in the exponential coupling models. 

In conclusion, we have performed a complete numerical study 
of interacting DE cosmologies for a few general types of time de- 
pendence of the DE-CDM coupling, concerning background evolu- 
tion, linear perturbations evolution, and nonlinear structure forma- 
tion. We have presented the first high-resolution hydrodynamical 
N-body simulations of structure formation in the context of vari- 
able coupling models to date. In our analysis, we have found that 
differently from the constant coupling case, halo density profiles 
and halo concentrations do not evolve in the same direction with 
respect to ACDM for all types of coupling evolution. In particular, 
depending on the type of background evolution determined by the 
coupling function, density profiles can be less overdense and cor- 
respondingly less concentrated than in ACDM, or vice versa. Fur- 
thermore, the growth of structures at large scales is also affected 
in a significantly different way according to the different types of 
coupling evolution. Finally, we find that the decrease of the halo 
baryon fraction already found for constant coupling models can be 
significantly enhanced in variable coupling cosmologies. Some of 
these effects alleviate tensions between astrophysical observations 
and the ACDM cosmology at small scales, and arise in cosmolog- 
ical models that contrarily to the constant coupling scenarios are 
not in stark conflict with present observational constraints on the 
background evolution of the Universe even in the presence of a sig- 
nificant coupling strength at low redshifts. Therefore, cosmologi- 
cal models with time dependent couplings in the dark sector might 
represent - for some specific forms of coupling evolution - a viable 
alternative to the standard ACDM concordance model. 
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